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PART  I.      PLAN  OF  RESEARCH 
I.     INTRODUCTION 

The  rational  design  of  pavements  must  lead  to  the  prediction  of  their 
performance  during  their  service  life.    The  overall  pavement  performance 
as  measured  by  its  serviceability  and  maintainability  under  induced  environ- 
mental and  loading  conditions,  in  turn,  depends  upon  the  structural  integrity 
of  the  component  layers.    The  engineering  properties  of  paving  materials  and 
geometrical  variables  such  as  thickness  and  relative  position  of  various  com- 
ponent layers  greatly  contribute  to  the  structural  integrity  of  pavement  systems, 
The  pavement  material  and  its  structural  arrangement  need  to  be  selected  to 
provide  optimum  serviceability  under  induced  load  and  environmental  factors. 
Recognizing  that  pavement  failures  can  be  placed  into  three  broad  categories 
of  durability,  stability,  and  fatigue  and  fracture,  bituminous  concrete  needs 
to  be  selected  and  designed  so  as  to  resist  detrimental  forces  of  load  and  en- 
vironment and  provide  satisfactory  performance. 

This  study  is  concerned  with  the  evaluation  of  fatigue  and  fracture  re- 
sistance of  bituminous  concrete.    It  is  aimed  at  an  investigation  of  the  char- 
acteristics and  behavior,  and  improvements  of  these  mixtures  as  related  to 
resistance  to  cracking.    As  it  has  been  shown  in  the  literature  review,  two 
approaches  to  such  a  problem  study  are  possible;  namely,  one  leading  to  a 
phenomenological  characterization  and  evaluation  of  asphaltic  concrete  mate- 
rial; and  the  other  based  on  the  mechanistic  approach  developed  at  Ohio  State 
University.    The  mechanistic  approach  utilizes  the  fracture  mechanics  concept 
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to  study  the  damage  processes  and  to  evaluate  the  fatigue  and  fracture  resist- 
ance of  asphaltic  mixtures. 

The  phenomenological  characterization  of  fatigue  in  bituminous  mixtures 
has  been  extensively  studied  within  the  past  few  decades.    In  a  recent  report, 
Saraf  (49)  has  presented  a  detailed  review  of  various  testing  and  analytical  pro- 
cedures used  for  the  phenomenological  characterization  of  fatigue  life  of  bitum- 
inous mixtures.    This  review  of  literature  indicates  that  various  methods  and 
analytical  procedures  have  been  employed  in  the  past  to  study  the  fatigue 
response  of  asphaltic  mixes  in  the  laboratory.    The  phenomenological  approach 
to  fatigue  is  credited  to  Monismith  (29-36)  s  Pell  (41-43)  and  others  (5,7,18)  in  which 
the  fatigue  life  of  a  flexible  pavement,  Nf,  is  related  to  the  maximum  tensile 
stress,  o",  or  tensile  strain,  e  ,  developed  in  the  under-side  of  the  bituminous 
layer  by  semi-empirical  relations  of  the  form: 


Nf   = 

oi(-) 

m9 
1         * 
C2(  7  ) 

for  controlled  stress  tests 

(1.1) 

Nf   = 

for  controlled  strain  tests 

(1.2) 

and 

where  c      c2,  m1 ,  m2  ,  are  constants  to  be  determined  experimentally  on 
simply-supported  or  cantilever  beams  using  prescribed  testing  procedures. 

According  to  these  equations,  for  a  repetitive  loading  of  a  controlled- 
stress  nature,  specimens  with  high  initial  moduli  or  stiffness  tend  to  perform 


*Underlined  numbers  in  parentheses  refer  to  the  list  of  references. 


most  satisfactory.    The  reverse  is  true  for  controlled- strain  fatigue  experi- 
mentation.   Therefore,  the  interpretation  of  the  test  results  and  the  selection 
of  mix  characteristics  required  for  pavement  design  are  dependent  upon  the 
mode  of  loading.    It  has  been  argued  (47)  that  since  the  constants  c  and  m  vary 
with  the  type  of  test  and  boundary  conditions,  they  cannot  be  considered  as  true 
material  constants. 

Despite  such  inherent  limitations,  the  phenomenological  approach  pro- 
vides a  reasonably  simple  procedure  which  has  been  universally  adopted  by 
various  research  organizations.    In  addition  to  the  variations  in  the  testing 
procedures  such  as  recommended  by  Monismith  et.al. ,  phenomenological 
evaluation  of  bituminous  concrete  has  also  been  carried  out  using  different 
geometrical  setups. 

Bazin  (2),  Savin  (50)  and  Coffman  (5)  have  used  trapezoidal-shaped 
cantilever  beam  samples  in  their  laboratory  investigations.    Jiminez  and 
Gallaway  (18)  used  an  apparatus  designated  as  a  "deflectometer"  and  a  rotating 
bending  fatigue  machine  has  been  used  by  Pell  (41-43)  in  this  study  of  asphaltic 
concrete  fatigue. 

Similarly,  researchers  have  argued  that  for  a  better  simulation  of 
fatigue  response  of  pavements  in  service,  the  mode  of  loading  needs  to  be 
selected  either  as  controlled  stress  or  controlled  strain  depending  upon  the 
pavement  thickness  and  other  geometrical  considerations  (1).    In  controlled 


stress  or  controlled  load  testing  as  represented  by  equation  1.1,  the  nominal 
stress  or  load  is  maintained  constant  throughout  the  duration  of  the  test.    If 
the  nominal  stress  level  is  maintained  constant,  the  testing  is  under  the  con- 
trolled strain  or  controlled  deflection  mode  (see  equation  1.2). 

The  available  experimental  results  have  indicated  that  the  fatigue  life 
of  a  given  sample  in  controlled  strain  tests  is  usually  higher  than  in  controlled 
stress  tests  (28,  29,_30)  when  compared  with  the  same  initial  nominal  stress. 

Most  of  the  previous  phenomenological  investigations  of  fatigue  response 
of  asphaltic  mixes  have  been  carried  out  using  simple  loading  histories,  in  which 
the  maximum  and  minimum  amplitudes  of  load  or  strain  were  constant  for  each 
cycle.    Although  this  criterion  has  provided  basic  information  regarding  the 
mechanism  of  fatigue  failure  in  a  given  sample,  it  is  far  from  the  true  loads 
sustained  by  pavements  in  service.    Compound  loadings  have  also  been  attempted 
by  Deacon  and  Monismith  (  7  )  to  study  their  effects  on  fatigue  response. 

Frequency  of  load  application  in  fatigue  tests  has  also  been  considered 
as  an  important  testing  variable.    Monismith,  et.al.  (27)  reported  that  the 
frequency  of  load  application  in  the  range  of  3.0  to  30.0  cycles  per  minute 
had  no  effect  on  the  specimen's  fatigue  behavior.    The  subsequent  work  of 
Deacon  and  Monismith  (  7  )  has  shown  that  the  increase  in  the  rate  of  loading 
significantly  decreased  the  fracture  life  for  the  type  of  test  employed  at  rates 
between  30  and  100  applications  per  minute. 


Raithby  and  Sterling  (42,43)  have  similarly  shown  that  the  rest  periods 
between  successive  loading  cycles  have  a  beneficial  effect  on  fatigue  performance, 
both  by  increasing  the  resistance  to  cracking  and  by  reducing  the  rate  of  loss  of 
dynamic  stiffness  due  to  repeated  loading.    Rest  periods  on  the  order  of  one 
second  increased  the  number  of  cycles  to  failure  by  a  factor  of  upto  5,  when 
compared  with  the  life  under  continuous  sinusoidal  cyclic  loading.    The  improve- 
ment in  life  was  less  at  high  temperatures;  it  also  appeared  to  be  influenced  by 
the  magnitude  of  the  applied  cyclic  stress,  although  this  effect  was  not  clearly 
established.    A  comparison  of  fatigue  performance  under  square,  sinusoidal, 
and  triangular  wave  forms  indicated  some  significant  differences,  but  these 
were  small  compared  to  the  effects  of  rest  periods.    It  should  be  noted  that  the 
effect  of  rest  periods  on  the  crack  growth  process  has  been  recently  investigated 
by  Majidzadeh  and  Kauffmann  (26)  by  testing  beams  on  elastic  founcations  with 
rest  periods  of  0.0,  0.4  and  0.8  seconds.    They  reported  that  for  such  test 
conditions,  there  were  no  significant  effects  of  rest  periods  on  the  crack 
growth  process. 

In  summary,  it  should  be  pointed  out  that  although  the  phenomenological 
approach  has  gained  universal  acceptance,  it  bears  the  limitation  that  it  cannot 
take  into  account  the  crack  initiation  and  propagation.    Nor  can  it  fulfill  the 
mechanistic  model  representation  of  fatigue  and  fracture  of  bituminous  or 
other  paving  materials. 


To  take  into  account  the  role  played  by  crack  propagation  and  the 
consequent  redistribution  of  stresses  within  a  layered  pavement  system,  a 
mechanistic  model  was  developed  at  Ohio  State  University  (23), (27),  (47)  in 
the  late  1960 's.    It  was  considered  essential  that  the  design  strategy  should 
incorporate  the  formation  and  development  of  failure  mechanisms  and  describe 
such  processes  involved  in  terms  of  invariant  material  properties,  loading, 
geometry  and  boundary  conditions. 

The  mechanistic  scheme  is  a  rational  strategy  which  utilizes  the 
principles  of  fracture  mechanics  to  explain  the  mechanism  of  damage,  the 
progression  of  crack  growth,  and  the  prediction  of  the  fatigue  life  of  pavements. 
Such  a  design  scheme,  as  shown  in  Figure  1,  is  based  on  the  postulate  that  the 
fatigue  life,  Nf,  can  be  described  by  the  process  of  crack  initiation,  growth 
and  ultimate  fracture.    The  material  parameters  required  for  a  detailed 
fatigue  and  fracture  characterization  are  parameters  associated  with  and 
identified  by  these  three  damage  processes: 

Nf   =    F  [c0,  A,  n,  K,  KIC]  (1.3) 

The  parameter  Kj£  is  the  fracture  toughness  describing  the  final  stages  of  crack 
propagation.    The  crack  growth  process,  however,  has  been  shown  to  be  ex- 
pressed in  terms  of  A,  n  and  K  by  either  a  simplified  power  law  equation: 

dc/dN    =    A  Kn  (1.4) 
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or  in  a  generalized  form  of 


ill  ^2  n3 

dc/dN   =    A1K      +  A2  K      +  A3  K  (1.5) 


where  Aj  ,  nj    =    material  constants  and  K  =  stress-intensity  factor.    The 
critical  value  of  the  stress-intensity  factor,  Kjp>,  is  a  material  constant  which 
describes  the  failure  criterion  for  both  fatigue  and  fracture. 

The  analysis  of  fatigue  and  fracture  of  pavement  systems,  using  the 
mechanistic  method,  as  shown  in  Figure  1,  is  a  two-step  process;  that  is, 
(i)  the  material  characterization  to  obtain  parameters  describing  the  crack 
resistance  of  paving  mixtures  and  (ii)  solution  of  boundary- value  problems 
to  obtain  an  estimation  of  the  stress-strain  distribution  within  the  layered 
system.    The  previous  research  at  the  Ohio  State  University  (26)  and  the 
theoretical  development  of  fracture  mechanics  at  other  institutions  have  pro- 
vided various  methods  of  solution  for  crack  problems.    The  computer  program 
and  analysis  procedures  presented  in  research  reports  (26)  (27)  are  typical  of 
such  analysis  procedures  available  for  the  design  of  pavement  systems. 

With  respect  to  material  characterization  task,  the  proposed  mech- 
anistic analysis  procedure  requires  the  determination  of  material  constants 
affecting  fatigue  and  fracture.    It  also  requires  the  sensitivity  analysis  of 
these  material  constants  as  affected  by  pertinent  testing  variables. 

As  an  example,  the  previous  studies  have  indicated  that  the  prediction 
of  the  fatigue  life  of  pavements  from  laboratory  tests  is  independent  of  the 


mode  of  loading  or  specimen  geometry.    Failure  by  cracking  is  defined  in  terms 
of  the  stress-intensity  factor  and  catastrophic  failure  by  the  critical  stress- 
intensity  factor.    The  results  of  these  studies  carried  out  at  the  Ohio  State 
University  have  also  indicated  that  the  crack  propagation  law,  in  a  more 
formalized  form,  is  applicable  to  asphaltic  mixtures  tested  in  various  geo- 
metries, boundary  conditions,  temperatures  and  modes  of  loading.    The  effect 
of  random  block  loading  and  the  levels  of  variable  amplitude  on  the  crack  pro- 
pagation has  also  been  formulated  (26).    The  mathematical  models  incorporating 
the  effect  of  plastic  yielding  at  crack  tip  locations  have  also  been  investigated. 
The  results  of  these  researches  have  led  to  the  following  general  conclusions: 

(i)    The  crack  propagation  model  dc/dN  =  A  Kn  satisfactorily  explains 
the  fatigue  performance  of  bituminous  material,  the  stress-intensity  factor 
being  the  dominant  parameter  controlling  the  crack  growth. 

(ii)    Theoretical  solutions  have  been  developed  that  accurately  predict 
the  stress  distribution  in  cracked  bodies  and  the  stress-intensity  factor  for 
beams  on  elastic  foundations,  slabs  on  elastic  foundation,  etc.    The  agree- 
ment of  theoretical  stress-intensity  factors  and  the  experimental  values  have 
been  found  to  be  excellent. 

(iii)    The  fatigue  crack  propagation  process  in  bituminous  materials 
is  considered  as  a  stress  interaction  type  phenomenon. 

(iv)    For  the  range  of  variables  investigated,  it  was  found  that  the  rest 
period  has  little  or  no  effect  on  the  fatigue  life  of  pavements. 


(y)    The  sequence  of  load  application  has  a  significant  effect  on  the 
fatigue  life  of  asphaltic  materials,  causing  a  delay  in  the  rate  of  crack  pro- 
pagation. 

(vi)    The  fatigue  life  of  pavement  slabs  subjected  to  loads  of  variable 
amplitude  can  be  predicted  from  tests  on  beams  resting  on  an  elastic  solid. 

(vii)    The  constant  n  in  the  power  law  dc/dN  =  A  K    varies  depending 
upon  the  material  characteristics  such  as  asphalt  content,  gradation,  etc. 

(viii)    The  mechanistic  approach  to  predict  cracking  of  flexible  pave- 
ments is  applicable  in  the  range  of  temperatures  from  41°  to  90°F. 

(ix)    The  critical  stress-intensity  factor,  KTC,  is  the  failure  criterion 
for  low  temperatures. 

(x)    The  concepts  of  linear  fracture  mechanics  can  be  used  to  predict 
the  fatigue  life  of  pavements  by  means  of  the  analysis  encompassed  by  the 
mechanistic  approach.    The  application  of  this  method  to  a  typical  flexible 
pavement  has  been  demonstrated  in  previous  studies  (27_,28) . 

To  utilize  this  mechanistic  model  for  pavement  design  purposes,  it 
is  required  that  the  material  characteristics  responsible  for  the  fracture 
and  fatigue  resistance  of  asphaltic  mixtures  to  be  determined  in  the  labora- 
tory under  field  simulated  conditions.    In  this  respect,  therefore,  this  study 
has  been  aimed  at  the  application  of  mechanistic  concepts  to  optimization  of  mix 
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variables  and  improve  the  fatigue  and  fracture  life  of  expectancy  of  the  bituminous 

mixtures.    To  achieve  such  objectives,  it  should  be  demonstrated  that  fracture 

and  fatigue  principles  are  applicable  to  a  variety  of  bituminous  mixtures  and  can 

provide  an  estimate  of  the  crack  resistance  characteristics  of  these  materials. 

In  addition,  the  analytical  and  experimental  fracture  mechanics  to  be  used  for 

mixture  optimization  must  be  simple  enough  to  lend  itself  to  routine  material 

characterization    and  mix  design  purposes.    It  must  also  recognize  the  effects 

of  variability  in  the  material  characteristics,  such  as  aggregate  size  distribution, 

asphalt  content,  and  other  mixture  variables.    The  scope  of  this  study  includes 

the  following. 

II.    SCOPE 

A.    Selection  of  Testing  and  Analytical  Procedures  for  Fracture  Resistance 

This  phase  of  the  study  involves  a  detailed  review  of  literature  and  pre- 
paration of  a  state  of  the  art  report  of  the  applicability  of  fracture  mechanics 
concepts  to  the  evaluation  of  fatigue  and  fracture  resistance  characteristics 
of  bituminous  concrete  materials. 

It  is  also  aimed  at  selecting  the  most  promising  existing  laboratory  test- 
ing procedures,  modified  as  necessary,  and  to  verify  these  procedures  through 
laboratory  testing  of  bituminous  concrete  mixtures.    Consideration  shall  be 
given  to  the  following  items  during  verification:    specimen  type,  size,  geo- 
metry, and  method  of  support;  mixture  type  and  maximum  size  of  aggregate; 
and  test  temperature. 
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With  respect  to  the  testing  temperature,  the  study  will  investigate  the 
agreement  between  elastic  and  viscoelastic  analysis  at  77  °F,  to  measure  the 
effect  of  time  dependent  material  response  on  fatigue  crack  propagation. 

To  evaluate  the  results  and  test  procedures  best  suited  for  crack 
resistance  studies,  analytical  techniques  such  as  finite  element,  boundary 
collocation  or  other  new  and  innovative  techniques  will  be  used  to  determine 
which  test  procedure  best  defines  crack  resistance. 

B.  Sensitivity  of  Various  Mix  Parameters  and  Evaluation  of  Formulations 
for  Improvement  of  Fracture  Resistance 

This  phase  of  the  study  involves  the  evaluation  of  the  effect  of  asphalt 
source,  consistency  and  amount,  aggregate  grading,  type  and  maximum  size; 
air-void  content  and  type,  and  the  amount  of  conventional  fillers  on  the  fracture- 
resistant  qualities  of  the  mixtures. 

This  phase  -will  also  include  an  evaluation  of  the  effects  of  polymeric 
additives  to  the  asphalt  and  the  effects  of  admixtures  such  as  asbestos,  rubber, 
sulfur  and  synthetic  constituents  on  fracture  resistance. 

C.  Recommendations  for  Routine  Mixture  Design  and  Testing  and  Suggested 
Formulations  for  Fracture-Resistant  Mixtures 

It  is  the  aim  of  the  study  to  recommend  a  mixture  design  procedure 

that  will  enable  optimization  of  fracture-resistant  qualities  of  typical  mixtures 

which  can  be  used  on  a  routine  basis.    Mixtures  formulations  required  for 
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optimizing  the  fracture-resistant  characteristics  of  typical  mixtures  will  also 
be  recommended  using  the  results  of  the  sensitivity  analysis  of  mixtures  vari- 
ables as  developed  in  Phase  n. 
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PART  n.      A  STATE-OF-THE-ART  REVIEW 
I.    FRACTURE    MECHANICS 

A.  A  General  Review 

One  of  the  important  aspects  of  modern  structural  analysis  is  the 
selection  of  a  criterion  of  failure  which  could  provide  an  accurate  estimate 
of  factors  affecting  failure.    The  traditional  design  criterion  often  does  not 
account  for  the  existence  of  flaws  and  inherent  defects  causing  premature 
structural  distress.    The  classical  approach,  which  is  based  on  the  selection 
of  a  limiting  applied  stress  as  compared  with  the  material  yield  stress,  <r 
might  not  always  guarantee  a  fail-safe  design  strategy. 

It  has  been  shown  that  the  theoretical  strength  of  materials  which  de- 
pends on  the  forces  of  molecular  cohesion  are  many-fold  larger  than  actual 
observed  values  of  strength.    It  is  obvious  that  material  defects,  imperfec- 
tions within  the  crystalline  structures,  and  other  forms  of  flaws  resulting 
from  manufacturing  and  handling  are  responsible  for  such  a  discrepancy.    The 
existence  of  such  flaws  as  joints,  cracks,  etc. ,  cause  a  redistribution  of 
stresses  and  stress  concentration  in  the  vicinity  of  structural  discontinuities. 
The  high  elevation  of  stresses  at  such  localized  regions  of  the  body  can  often 
result  in  catastrophic  failure  of  the  structure,  even  at  normal  stress  levels 
much  less  than  the  yield  strength  of  the  material. 
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To  account  for  such  discrepancies  between  the  observed  and  theoretical 
failure  limits,  a  classical  theory  of  fracture  mechanics  has  been  set  forth  to 
provide  a  rational  and  refined  analysis  of  the  degradation  of  the  strength  due 
to  inherent  flaws.    In  brief,  in  the  fracture  mechanics  approach,  the  existence 
of  flaws  are  assumed  to  be  responsible  for  the  elevation  of  stresses  at  crack 
tip  locations.    When  the  applied  load  is  increased,  the  stresses  around  the 
flaws  reach  a  limiting  value  and  fracture  becomes  possible. 

The  basis  of  linear  fracture  mechanics  is  the  paper  by  Griffith  (13)  which 
was  unrecognized  for  thirty  years.    In  1920,  Griffith  proposed  that  brittle  bodies 
fail  because  of  the  presence  of  numerous  internal  cracks  or  flaws  which  pro- 
duce local  stress  concentrations.    He  also  stated  that  the  elastic  body  under 
stress  must  transfer  from  an  unbroken  to  a  broken  state  by,  a  process  during 
which  a  decrease  in  potential  (elastic)  energy  takes  place.    He  postulated  that 
fracture  instability  is  reached  when  the  increase  in  free  surface  energy  (surface 
tension),  caused  by  the  extension  of  the  crack,  is  balanced  by  the  release  of 
elastic-strain  energy  in  the  volume  surrounding  the  crack.    The  Griffith  equa- 
tion for  propagation  of  an  internal  crack  of  length  2c    in  an  infinite  thin  plate 
under  a  uniform  stress,  <r   ,  normal  to  the  plane  of  the  crack,  is: 

*2  >-  *&- 

where  y   =    the  specific  surface  energy 

E    =    the  modulus  of  elasticity 
c    =   crack  length. 
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In  many  engineering  materials  such  as  ceramics  and  glass,  fracture 
occurs  with  hardly  any  deformation.    However,  in  metals  and  other  engineer- 
ing materials,  plastic  deformation  always  takes  place,  and  Griffith's  classical 
surface  tension  concept  is  not  suitable.    Both  Orowan  (52,53)  and  Irwin  (17, 18, 19) 
independently  came  to  the  conclusion  that  the  slight  plastic  flow  which  occurs 
in  the  brittle-fracture  case,  absorbs  a  large  amount  of  additional  energy  re- 
quired to  create  new  surfaces. 

Mathematically,  Griffith's  criteria  for  crack  propagation  is: 

SU       >       5UST  (2.2) 

where  5U  is  the  decrease  in  potential  energy  due  to  increased 

crack  surface, 
5UgT       is  the  increase  in  surface  energy  due  to  increased 
crack  surface. 

Irwin's  criteria  for  crack  propagation  can  be  represented  as: 


5U      >     5UgT    +         5UpL  (2.3) 


where  SUpi,      is  the  plastic  energy  dissipated  due  to  increased 

crack  surface. 

Irwin  recognized  that  the  plastic  energy  dissipated  is  much  larger 
than  the  surface  energy  dissipated  and  therefore,  proposed  to  ignore  the 
latter. 
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The  terms  above  involving  increase  in  energy  quantities  are  evaluated 
with  respect  to  the  increased  crack  surface   5  A    and  hence,  the  change  in 
energies  could  be  referred  to  as  rates  of  energy  input  or  dissipation. 

The  theory  by  Griffith  on  energy  balance  and  its  subsequent  modification 
by  Irwin  and  Orowan  are  necessary  conditions  for  the  onset  of  fracture.    As 
such,  this  theory  cannot  conveniently  characterize  all  types  of  fracture  ob- 
served in  fracture  tests.    Irwin  thus  proposed  that  the  local  stress  field  sur- 
rounding the  crack  tip  be  used  in  place  of  the  total  input-output  energy  rate 
for  such  characterizations. 

The  relationship  between     5U     ,  the  input  energy  rate,  and  the  local 

stress  field  can  be  obtained  by  considering  the  situation  where  a  short  segment , 

8x,  of  a  two-dimensional  crack  is  closed  by  imposing  a  force,  <r      (-  8x,0), 

on  the  crack  surface  as  shown  in  Figure  2.1  (21) (32).    The  total  strain-energy 

absorption  rate  in  this  reverse-loading  problem  is  equal  to    8  U^.  ,  which  in 

turn  is  equal  to  the  strain-energy  release  rate,  G,  for  a  crack  extension  of 

8  x,  or 

5x 

81V       =     G    •      Sx      =     J      Uy(0,0)    •    <r*y(-  5x,0)  dx        (2.4) 


where  uy(0,  0)  is  the  crack  closing  displacement  when  the  crack  closes  from  the 


origin  of  the  x-y  coordinate  for  a  length  of  8  x. 
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<ry;(-s*,o) 


Figure  2.1.    Closing  of  a  Crack  Tip 

It  can  be  assumed  further  that  the  plastically  yielded  region  ahead 
of  the  crack  tip  does  not  change  the  state  of  stress  significantly.      As  a  re- 
sult,    ^(0,0)  and    <r  *  (-    5x,0)  can  be  obtained  by  the  elastic  stresses  and 
elastic  displacements  obtained  through  the  known  elastic  state  in  the  vicinity 
of  the  crack  tip.    In  particular,    1^.(0,0)  is  the  crack-closing  displacement 
with  the  crack  tip  located  at  the  origin  of  the  coordinates  (0,0)  and  <r  £v(-  5x,0) 
is  the  associated  closing  stress  which  is  equal  to  the  normal  stress  ahead  of 
the  crack  tip  located  at  (-  §x,0).    Therefore,  all  problems  in  linear  fracture 
mechanics  can  now  be  converted  to  problems  in  linear  elasticity  for  which 
solutions  are  possible. 

It  is  also  shown  that  all  linear  elastic  solutions  to  problems  involving 
cracks  in  homogeneous,  isotropic  material  show  an  inverse  square  root  de- 
pendency of  the  crack  tip  singularity. 
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Considering  the  "Mode  I"  crack  tip  deformation  problem,  as  shown  in  Fig- 
ure (2.1),  Westergaard's  stress  function,  frequently  used  to  solve  two-dimensional 
problems  in  cracked  structures  in  the  vicinity  of  the  crack  tip,  is: 

zi   -      Vp  {2'5) 

where    p  =    re^    (in  polar  coordinates. 

As  p  >  0,  however,  the  complex  analytical  function,  f  (p),  approaches 
a  real  constant,  — j=  ,  where  Kj  is  the  opening  mode  of  stress-intensity  factor 
and  is  defined  as: 


Kj   =      lim       <ryy(r,  6  =  0\    •     V2  7Tr  (2.6) 

r— ►O 

The  local  stresses  in  terms  of  local  coordinates  for  Mode  I  fracture  are: 


XX 


.    e     .    se 

1  -  sm  ~2~   sin  ~2~ 

KI  0        J   ,         .      d      .     W 

yy  >     =        ~~ZZZ —     cos  ~2~      )        sin  2    sin  ~?~ 

V2  ?rr 


°"xy 


sin  -jr-     cos  -^~ 

Similarly,  the  local  displacement  components  in  the  x  and  y  directions, 
ux  and  Uy,  can  be  expressed  in  terms  of  Kj  as; 


uy 


cos-|-  (X-  1    -  2  sin2-|- 


i1  +v)     Kj  r2  (2.8) 

V2~7T     E  sin  A  (X+ 1  -  2  cos2-^-) 

2  2 


where  X   =  3  -  4p    for  plane  strain 

=  (3  -  V)/(l  +V)  for  plane  stress 

V  =  Poisson's  ratio 

E    =  the  modulus  of  elasticity. 

26 


The  substitution  of   $  -  7T  in  Eq.  (2.8)  gives  ux  =  0  and 

1 
Uy(0,0)    f    (1  -0  r2       •     Kj(A-+l)  (2.9) 

or  4  (1  -    V2)        i 

uy  (0,0)    =    \/G^T    F  r2    KI   f°r  plane  strain  (2.10a) 

and 

u    (0,0)    =  4  r2    Kj   for  plane  stress  (2.10b) 

V2lF   E 
* 
The  imposing  stress    cr       (-   5x,0)  is  obtained  by  shifting  the  origin  in 

Figure  (2.1)  from  (0,0)  to  (-  8x,0)  and  evaluating  cr      in  Eq.  (2.7b)  for  r  =    8x-  r, 

9  =  0.    This  gives,         *  kt 

o-      (-  8x,0)    =        — (2.11) 

yy  V2  IT    (  8  x  -  r) 

The  strain-energy  release  rate  due  to  a  crack  extension  of    8  x  in  a  state 
of  plane  strain  is  then  given  by 


2  (1  -  P  ) 


2s  „  8x 


8U-  =  G*  6x  =  ~7T  '  Ki      /      ^T      dr 


o  8x     / 

■  J  ( 


(2.12) 


which  upon  substitution  r  =     8  x  sin        and  evaluation  of  the  integral  gives 

r      *              (1  -  V2)      „2  . 

Gj   8x   = Kj     •      8x 

or  (1  -    V2)         o 

Gj  =    g Kj       (plane  strain)  (2.13) 

For  a  plane  stress  condition,  a  similar  relation  can  be  obtained  by  simply  substitute 
Eq.  (2.10b)  instead  of  Eq.  (2.10a)  in  Eq.  (2.4).    This  gives  Gj  =  Kj/E  for  plane  stress, 
The  subscript  I  denotes  the  opening  mode  fracture. 

At  the  onset  of  fracture,  the  elastic  energy  release  rate,  following  the 
Griffith-Irwin  theory,  can  be  replaced  by  an  equivalent  material  constant  desig- 
nated by  the  critical  strain-energy-release  rate  or  GTp.    Since  the  stress  in- 
tensity factor  is  directly  related  to   Gj  ,  the  Griffith-Irwin  theory  can  also 
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be  restated  in  terms  of  the  stress-intensity  factor,  Ky,  which  becomes  a 
material  constant  Kjq  at  the  onset  of  rapid  fracture.    This  critical  stress- 
intensity  factor,  KTr,,  is  referred  to  as  the  fracture  toughness  of  the  material 


and  has  the  dimensions  of     (stress     •      vTength  ) 

With  regard  to  a  general  state  of  fracture  in  which  in-plane  sliding 
(Mode  II)  and  tearing  (Mode  III)  are  present  (Figure  2.2),  the  stress  distri- 
bution at  the  end  of  cracks  can  be  written  as: 


_i 
cr,,     =     r  2 


Kifij(0)  +    Knf"(^)    +    Knifijn(0 


+     other  nonsingular  terms. 


(2.14) 


where  r,  0       are  polar  coordinates  introduced  at  the  crack  tip 

(Figure  2.3) 

Kj,  Kjj,  Kjjj  and  corresponding  symbols  I,  II,  and  III  reflect 
three  modes  of  fracture. 


As  was  pointed  out,  the  terminal  state  of  fracture  can  be  represented 
by  the  critical  value  of  stress-intensity  as  designated  by  Kjq.    The  interrela- 
tion of  Irwin's  stress  intensity  factor  K_  and  Griffith's  energy  release  rate  Gj 
are  represented  in  Equation  (2. 13).    Although  similar  relations  between  Gjj 
and  Ktt  ,  and  K  TT  and  Grjjcan  also  be  written,  the  physical  meanings  of  Kjjq 
and  Kjjjc  have  not  been  well-defined.    Therefore,  the  Mode  I  fractures  have 
remained  a  significant  method  of  analysis.    The  Mode  I  critical  stress-intensity 
factor,  Kj£,  as  designated  by  fracture  toughness,  has  also  been  considered  as 
an  important  engineering  tool  which  provides  the  designer  with  an  additional 

material  parameter  for  crack  resistance  evaluation. 
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The  difference  between  Kj  and  Ky^  terms  is  quite  significant  and  should 
be  clearly  understood.    The  term  Kj,  or  a  stress-intensity  factor,  is  determined 
analytically  and  is  a  function  of  load,  geometry  and  crack  size.    Hence,  the 
analytical  form  and  the  magnitude  of  KT  varies  from  one  system  to  another, 
depending  on  load  and  geometrical  variables.    On  the  other  hand,  the  fracture 
toughness,  Kj£,  is  a  material  constant,  independent  of  crack  geometry,  loading 
conditions  and  other  physical  variables  and  is  analogous  to  strength. 

It  should  be  noted  that  the  K~,  concept  is  primarily  restricted  to  labora- 
tory  and  material  characterizations  where  the  Mode  I  fracture  is  applicable. 
That  is,  the  Griffith-Irwin  theory  is,  in  fact,  a  scalar  theory,  in  that  only 
the  critical  values  of  a  scalar  Gjq  and/or  KjC  are  known  (73  ).     The  direc- 
tion of  crack  propagation  is  assumed  to  be  normal  to  the  load  and  the  crack 
front  must  be  straight  (Figure  2.4). 

In  reality,  however,  in  most  structural  components  the  flaws  and  cracks 
are  seldom  aligned  perpendicular  to  the  direction  of  load.    Such  a  deviation  in- 
validates the  classical  Kc  theory  of  fracture,  in  which  cracks  must  always 
be  normal  to  applied  tensile  stress.    As  a  departure  from  this  classical  theory, 
Sih  (73)  has  presented  a  theory  of  fracture  based  on  the  field  strength  of  the 
local  strain-energy-density.     This  theory  has  the  inherent  advantage  of  being 
capable  of  treating  all  mixed  mode  fracture.    The  mixed  mode  fracture  analy- 
sis is  of  particular  interest  in  the  study  of  bituminous  materials.     It  has  been 
shown  that  for  mixtures  with  large  size  aggregate,  when  tested  as  a  beam  on 

elastic  foundation  or  in  real  pavement  systems,  the  crack  surface  follows  a 
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zig-zag  pattern  and  as  a  result,  the  stress  component     cr        responsible  for 
the  crack  growth  would  no  longer  remain  perpendicular  to  the  crack  surface. 

Under  such  conditions,  the  stress  state  ahead  of  the  crack  will  be 
governed  by  at  least  two  stress-intensity  factors  (K    and  K   )  instead  of  a  simple 
parameter,  IC.      The  elastic  stress  distribution  at  the  crack  tip  can  be  written 
as: 
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where 


(2.15) 


K, 


KjA/SF 


K2    =     K^ 


Under  such  a  state  of  stress,  the  failure  condition  can  be  evaluated 
using  the  energy  concept.    Considering  that  strain  energy  of  the  system  in  the 
crack  vicinity,   <p  ,  can  be  expressed  in  terms  of  stresses,  for  an  element 
r  A   d     ,   Ar  ahead  of  the  crack  tip,  the  strain  energy  function  is  (Figure  2.  3). 
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Then,  by  inserting  the  stress  components  from  the  previous  equation  (2.15), 


d<f>     =     1 
dA  r 


A     K  2    +    A     K  2    +    2A 


11    1 


22    2 


12     1       2 


(2.17) 


where  the  term  in  the  bracket  is  defined  as  strain-energy-density  factor  S, 


S     =     AnKx     +    A22K2      +    2A12K1K2 


(2.18) 


A..  1 ,  A1?  and  A22  are  constants  related  to  the  elastic  parameters  and  element 
direction  as  given  by: 
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IE 


4(1  -  v  )(i  -  cos  $  )  +  (1  +  cos  0)(3  cos0   -  1) 


For  the  condition  that  all  three  modes  of  fracture  are  present,  S  is  written  as: 


S     =     AUKX2     +     A22K2      +     A33K32     +     2A19K1K9 


l12^i^2 


where 
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If  the  loading  condition  and  crack  propagation  is  such  that  Mode  I 
analysis  applies,  then: 
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or 


(1  +  y)  (1  -  2v)  2 

S    =  2  77E  KI 


(2.21) 


When  unstable  crack  extension  or  failure  occurs,  the  parameter    S 
approaches  its  critical  value,  given  by  Sc: 

(1  +  v  )(1  -  2v)  2 


2  n  E 


K 


IC 


(2.22) 


where  in  this  case,  S    is  a  material  constant. 


Under  an  in-plane  shear  condition,  where  only  the  Mode  II  fracture  is 
possible,  the  local  stress  field  contains  K^  (K£  =  ^ii/\/1t~ )  alone,  and  the 
strain-energy-density  is  written  as: 


S     = 


1  +v 

8TTE 


4(1  -  v  )(1  -  cos  6  )    +    (1  +  cos  6  )(3  cos  0-1) 


K. 


II 


where 


K 


II 


t  y/Wc 


(2.23) 


in  which     r   is  applied  shear  stress  and   c    is  the  crack  length.    At  the  point 
of  crack  instability,  Sih  has  shown  that 


6E 
6  IT  E 


2  (1 


-„)-     -2] 


^e 


2  (1  -v)  -    v2\   •     Kj2IC 


(2.24) 
(2.25) 
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Since  Sc  is  a  material  constant,  it  can  be  argued  that  Equation  (2.24^  and 

Equation  (2.25)  must  be  equal;  and  thus  setting  these  equations  equal  to  each 

2 
other  and  solving  for  Kjr,-,    , 


K 


2  3  (1  -  2  v)  2 

*      KIC  (2.26) 


IIC 


2  (1  -  v)  -  v2 


A  similar  relation  can  be  developed  between  Kjq  and  Kjjtq  ,  when  considering 
the  Mode  III  fracture  process. 

The  preceding    relations  are  derived  from  Sih's  multi-model  frac- 
tum  concept  which  is  based  upon  the  postulation     that  crack  spread  occurs 
in  the  direction  of  maximum  potential  energy  density  of  the  system.      Con- 
sidering the  potential  energy  is  equal  to  the  negative  of  the  strain  energy  S, 
therefore  the  necessary  and  sufficient  condition  for  crack  growth  is  the  min- 
imum of  the  strain  density, S.    The  angle  0  =   8Q,  which  makes  S  a  minimum, 
determines  the  angle  at  which  the  crack  propagates.      That  is,  the  crack  ini- 
tiation occurs  in  a  direction  determined  by  the  stationary  value  of: 


a  S 


0    when    0      =      dQ  (2.27) 


The  strain  energy  density  factor,  S,  can  be  resolved  into  two  components, 
one  responsible  for  the  change  of  shape  as  S^  and/or  the  change  in  volume 
S   ,  that  is: 

S      =       S,     +     S  (2.29) 

a  v 
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The  crack  initiation  occurs  under  conditions  of  Smjn  for  which  Sih  (73)  has 
shown  that  fracture  occurs  along  the  plane  where    Sv   >  S^.      On  the  other 
hand,  the  direction  6q  =  cos"     (1  -2v)  along  which  Smax  occurs  (Mode  I  only), 
corresponds  to  maximum  yielding  on  the  large  dimension  of  the  plastic  zone. 

For  the  analysis  of  fracture  problems  involving  large  plastic  deforma- 
tions, the  deformation  theory  of  plasticity  has  been  used  by  Rice  (66),(67)  to 
arrive  at  a  parameter  needed  for  characterization  of  crack  tip  area.     The 
path  independent  J  integral  is  a  parameter,  the  value  of  which  depends  on  the 
near  tip  stress-strain  field.     The  path  independency  of  J  integral  has  been 
shown  to  be  valid  for  linear,  non- linear  elastic,  and  plastic  materials  under 
monotonic  loading  condition.    The  basis  of  this  formulation  is  that  the  path- 
independent  nature  of  the  integral  allows  integration  to  be  carried  out  along 
the  path  away  from  the  crack  tip  instead  of  along  a  region  close  to  the  tip. 

The  J  integral  is  defined  as: 


/ 


W  dy    -     T    (|^-  )  d  s  (2 .  30) 


r 

where    T    is  any  contour  surrounding  a  crack  tip,  W  is  strain  energy  density 


(Fig.  2.5): 


/emn 
°i]  d  eij 


W      =  [2.31 

0 

u,  displacement  vector,    T,  traction  vector,  defined  by  the  outward  normal 
to  the  contour.      Rice  (66)  has  pointed  out  that  the  crack  tip  deformation  and 
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Figure  2.5.    Crack  Tip  Coordinate  System  and 
Arbitrary  Line  Integral  Contour 
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energy  conditions  are  reflected  by  the  J  integral.    The  J  integral  might  be  inter- 
preted as  the  potential  energy  difference  between  two  identically  loaded  bodies 

having  neighboring  crack  size  c  and  c  +  Ac.    That  is, 

J  =  -  dU 


dc  (2.32) 

where  U  is  the  potential  energy,  and  c  is  the  crack  length.    In  analogy  with  the 
linear  crack  fracture  mechanics,  the  area  between  two  monotonic  load-deflection 
curves  (Fig.  2.6)  for  cracks  of  c  and  c  +  Ac  is  defined  as  Area  =  OABO.    That  is, 
J  can  be  determined  from  experimental  load-deflection  curves  of  varying  initial 
crack  length.    It  should  be  noted  that  since  the  deformation  theory  of  plasticity 
does  not  permit  unloading  conditions,  the  J  integral  is  only  applicable  to  crack 
initiation  rather  than  crack  propagation. 

For  a  linear  elastic  case,  J  is  identical  to  G,  the  energy  release  rate, 

and  may  be  expressed  as 

P2  dL  K2 

J   "    G   "    2B     '    be      =      E'  (2*33) 

where  P  =  applied  load 

B  =  width  of  the  material 

L  =  material  compliance 

c  =  crack  length 

E'  =  E  for  plane  stress  and  E/(l  -  f2)  for  plane  strain. 

For  a  rigid-plastic  material,  the  deformation  or  displacement  A  is  un- 
limited at  the  limit  load,  P  =  PL,  whereas  for  P  <  P-r  ,    A  =  0.    In  this  case  J 

is  given  by 

A  *PL 

J  =    b    ■  "5T-  <2-34) 

where  the  derivative    dPL/dc  can  be  evaluated  experimentally.    The  J 
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integral,  therefore, can  be  used  as  a  parameter  required  for  characterization 
of  crack  tip  problems  in  both  linear  and  non-linear  systems. 

B.     Crack  Tip  Stress  Analysis  For  Pavements 

From  the  viewpoint  of  design  and  analysis  of  pavement    layered  sys- 
tems, the  Mode  I  as  well  as  mixed  mode  fractures  have  been  found  to  be  of 
great  utilization.    The  Mode  I  fracture,  or  K-Theory,  is  primarily  applicable 
to  fatigue  and  fracture  analysis  of  laboratory  specimens  as  well  as  design  of 
flexible  layered  systems.    The  Multi-Modal  fracture,  on  the  other  hand,  can 
be  used  to  study  the  crack  propagation  in  composite  pavements,  continuously 
reinforced  concrete  and  development  of  overlay  design  strategies. 

The  opening  mode  (Mode  I )  or  bending  mode  fracture  has  been  used 
extensively  to  study  the  fatigue  and  fracture  characterization  of  asphaltic 
mixtures.    The  research  works  at  The  Ohio  State  University,  on  the 
mechanistic  modeling  of  fatigue  process  of  flexible  pavements,  have  been 
primarily  concerned  with  the  laboratory  analysis  of  Mode  I  fracture  involving 
the  use  of  asphaltic  beams  or  slabs  resting  on  elastic  foundation. 

The  basic  step  in  such  a  method  of  analysis  is  the  development  of  an 
interrelation  between  Irwin's  local  stress-intensity  factor,  K,  and  the  crack 
size,  c.    A  normalized  (K/P  -  c/d)  relation  can  then  be  used  to  investigate 
the  crack  propagation  process  of  laboratory  specimens  subjected  to  bending 
mode  fatigue  loading. 
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Experimentally,  the  K-c  relation  can  be  determined  for  different  types 
of  loading  and  for  any  geometrical  and  crack  patterns  by  a  simple  compliance 
measurement  technique.    According  to  this  procedure,  the  compliance  L, 
crack  length  c,  and  the  rate  of  change  of  compliance  with  crack  length  (d  L/  ^  c 
versus  c)  can  be  determined  experimentally  using  prenotched  specimens.    Then 
the  relation  of  the  stress-intensity  factor  with  the    <3L/dc    compliance  can  be 
represented  by: 


/ 


E  dL 

K/P    =V  2  0-  ~    pl)    '  dc  (2.35) 


where  P  =  applied  load 

E    =  Young's  modulus 

v   -  Poisson's  ratio 

L   =  compliance  or  inverse  slope  up  the  load/deflection  diagram 

c    =    crack  length 

K   =    opening  mode  stress-intensity  factor. 
Such  a  procedure  has  been  employed  by  Majidzadeh,  et  al. ,  (33)  (34)  (35), 

for  both  beam  on  elastic  foundation  and  slab  supported  on  elastic  foundation. 
In  Figures  (2.7)  and  (2.8),  the  normalized  compliance  L/L0  relations  for 
asphaltic  beams  on  elastic  foundation  and  compliance  L  for  pavement  slabs 
resting  on  elastic  support  are  presented.    The  comparisons  of  the  K/P  -  c  rela- 
tions developed  using  experimental  compliance  procedures  and  numerical  tech- 
niques, boundary  collocation  and  finite  elements  for  beam  on  elastic  foundation, 
are  presented  in  Figure  (2.9). 
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As  was  indicated,  the  K/P  -  c  relation  can  be  obtained  using  such  ana- 
lytical procedures  as  finite  element,  boundary  collocation,  etc.    The  details 
of  such  programs  have  been  presented  in  References  (34)  and  (35).    These 
analytical  procedures  include  a  method  of  analysis  for  slabs  resting  on  elas- 
tic foundations,  beams  on  elastic  foundations,  and  a  prismatic  finite  element 
program  suitable  for  an  overlay  structural  analysis.    A  typical  K/P-c 
relation  developed  for  slabs  resting  on  elastic  foundation  is  shown  in  Figure 
(2.  10).    hi  Part  HI  of  this  report,  further  applications  of  such  procedures 
are  discussed. 

A  number  of  investigators  (74, 77)  have  utilized  the  K/P  -  c  relations 
developed  for  simply  supported  beams.    The  most  well-recognized  form  of 
such  a  relation  is  from  the  Winne  and  Wundt  equation  as: 

KI      =      °n    (1-"2)h   f(cM  (2.36) 

where  <*n     =     nominal  bending  stress  at  the  root  of  the  notch,  psi 

6M  , 

=  -tttj vtr-    ,  where 

b(d  -  c)^     ' 

M   =   bending  moment  inpsi.    =  Pl/4 

P    =   applied  concentrated  load  at  midspan ,  pounds 

1     =   span,  inches 

v  =  Poisson's  ratio 

b  =  width  of  the  beam,  inches 

c  =  crack  depth  ,  inches 

d  =  depth  of  the  beam,  inches 

h  =  d  -  c 

f  (c/d)  =  function  of  the  geometry  given  in  Figure  2. 11. 
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o 

The  f(c/d)  can  also  be  written  as  f(c/d)  =    nc    (1  -  c/d)    .     For  deeply 
notched  beams  f(c/d)  has  a  limiting  value  of  0.52  (see  Fig.  2.11).    An  analytical 
formulation  of  K  -  c  relation  by  Srawley  (74)  leads  to  a  normalized  relation  be- 
tween  Y    and  c/d  as  shown  in  Figure  (2.12)  from  which  the  opening-mode  stress- 
intensity  factor  for  three-point  bending  beams  can  be  evaluated. 

Similarly,  Gross  et.al.  (14),  utilizing  boundary  collocation  and  Green's 
function  methods,  have  developed  a  numerical  function  for  stress-intensity  factor 
for  simply-supported,  three-point  bending  beams  with  centrally  located  crack  as: 

K     =    <r  %/c~     F(c/d^  (2.37) 

where  cr    =    -^r  ,  M  =  Pl/4  . 

B  d^ 

For  1/d  equal  to  4(see  Figure  2.13), 

F(c/d)    =    1.090  -  1.735  (c/d)  +  8.20  (c/d)2  (2.38) 

and  for  1/d  equal  to  8  see  Figure  2.13, 

F(c/d)    =    1.107  -  2.120  (c/d)  +  7.71  (c/d)2  -  13.55  (c/d)3        (2.39) 
+  14.25  (c/d)4 

A  fracture  toughness  K  calibration  curve  has  also  been  recommended 

by  ASTM  using  Brown  and  Srawley' s  analysis  (  5 ) ,  as  given: 

6M  c2 

Kt    =    Y 

1  Bd2 

where,  for  a  three-point  bending, 

Y    =    1.96  -  2.75  (c/d)  +  13.66  (c/d)2  -  23.98  (c/d)3  +  25.22  (c/d)4 

(2 .  40) 
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Figure  2. 12    Dimensionless  Stress- Intensity  Factor  for  a 
Simply  Supported  Beam  (from  Ref .  14) 
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Also  c,  d,  B  are  crack  length,  beam  depth  and  beam  width  respectively.    M  is 
the  bending  moment,  either  in  3-point  or  4-point  bending. 

The  K  calibration  in  terms  of  c/d  ratio  and  for  3-point  bending,  may 
also  be  written  as: 


V 

KI  Bd3/2 


2.9(c/d)1/2    -4.6(c/d)3/2    +    21.8(c/dh5/2 


-    37.6(c/d)7^2    +    38.7(c/*9'2 


(2.41) 


where  P„  is  the  applied  load  determined  according  to  the  procedure  given  on 
page  52,  and  1  is  the  span  of  the  3-point  bend  specimen. 

Among  other  important  developments  of  the  analytical  and  mechanistic 
aspects  of  the  rational  pavement  design  are  the  analysis  of  crack  plate  on 
elastic  foundation  (64> ,  and  the  prismatic  finite  element  with  multi-modal 
fracture  analysis  (3S) .    In  either  case,  the  multi-modal  fracture  analysis  is 
carried  out  using  finite  element  method  and  the  strain  energy  density  S  and 
factors  Kj  -  c  and  Krj;  -  c  relationships  are  calculated,  relating  the  crack 
growth  rate  dc/dN  to  the  stress-intensity  factors. 
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II.    FRACTURE  TOUGHNESS 

A.  General  Review 

It  was  shown  previously,  as  a  result  of  flaws  and  imperfections  in 
materials,  that  there  exist  local  stress  concentrations.    The  elevation  of 
stress  and  its  redistribution  at  crack  tip  positions  are  reflected  by  the  stress 
intensity,  K,  and  strain  energy  release  rate,  G.    Physically,  the  parameters 
Kr,  Kjj,  and  K™.,  related  to  the  general  state  of  local  stress  field,    are  inten- 
sities of  load  transmittal  through  the  crack  tip;  similarly,  the  GT,  G    ,  and 
G      are  strain  energy  release  rates  for  three  modes  of  fracture. 

The  interrelation  between  the  stress  intensity  factor,  K,  and  strain 
energy  release  rate,  G,  are  written  as: 

For  plane  stress:  K2    =     EG  (2.42) 

For  plane  strain:  K2    =        E  G  (2 .  43) 

In  engineering  analysis  of  fracture  problems,  it  is  preferred  to  work 
with  the  stress-intensity  factor  rather  than  the  strain  energy  release  rate. 
The  reason  is  that  the  stress  intensity  K,  due  to  the  superimposed  effects  of 
various  stress  fields,  are  linear  additions  of  corresponding  K  values,  as 
given  by: 
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Kx   =    Kj(l)  +  Kj(2)  +  Kj(3)  +  .   .   .  +  Kj(n) 

KII   =    Kll^  +  Kn(2)  +  KH<3)  +  •   •  •  +  Kn(n)  (2.44) 

Km   =    Km(l)  +  Km(2)  +  Km(3)  +  .   .   .+Km(n) 

On  the  other  hand,  in  strain  energy  terms,  the  superimposed  field  is  given  by: 

2 


Gj  =      [GjCIP    +    Gj(2)2    +  .  .  .  +  Gi(n)2  J 

Gn  =      [Gn(1)*  +  Gn(2)*  +  •  •  •  +  Gn<n)^  ]  2  <2-45> 

Gm  =     [Gm{l)h  +  Gm<2>*  +  •  •  •  +  Gm(n)*  ] 2 


Secondly,  in  the  analysis  of  fracture  mechanics  problems,  considera- 
tion should  be  given  to  the  definition  of  plane  stress  and  plane  strain  conditions. 
In  the  field  of  mechanics,  the  plane  stress  condition  is  defined  as: 

a     =   0  .        T  yz   -  rxz   =   0 

z  ' 

and  the  plane  strain  is  defined  as: 

•  az    =         "  (    °x   +    ay)»      Txz    =Tyz    =   0  ,      ez    =   0 

In  constrast  to  general  field  of  mechanics,  the  fracture  mechanics  defi- 
nition of  plane  strain  and  stress  conditions  are  more  restrictive  and  are  only 
applicable  to  crack  tip  locations.    In  such  areas,  the  plastic  zone  ahead  of  the 
crack  tip  is  constrained  from  elastic  deformation  by  the  presence  of  elastic 
material  along  the  crack  front.    If  the  plastic  zone  is  small  as  compared  to 
the  dimension  and  the  length  of  the  crack  front,  a    plane   strain  condition 
prevails.    The  size  of  the  plastic  zone  can  be  estimated  for  a  plane  strain 
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condition  as: 


2  n  a 


(1-21/) 


2         P 

kt  cosz  e 

2"     )  2 


(1  -  2  »/  )2  +  3  sin2  9 
2J 


+  KT  Kn  sin  9 


3  cos  6  - 


+  K 


II 


3  +  sin2  e    ((I  -  2  V  )2  -  9  cos2  € 


2/J 


+  3K 


In  the  crack  plane,  where  0=  0  and  for  a  two-dimensional  case: 


r      =  1 

y      

2rra 


2" 


KT2  (l-2i;)2  +  3^KII2   +4^  ) 


For  Mode  I  fracture,  this  equation  reduces  to: 

,2        /  v,  \  2 


=    (1-21/V 
y  2;r 


(f) 


for  plane  strain. 


nil 


(2.46) 


(2.47) 


(2.48) 


As  the  material  body,  with  pre-existing  flaws,  are  stressed  to  failure, 
the  stress  intensity  factor  K  and  strain  energy  release  rate  G  are  increased 
approaching  their  limiting  values.    These  performance  limiting  parameters, 
for  the  plane  strain  condition  are  designated  as  Kj^  and  Gjq  which  have  been 
shown  that  they  can  be  taken  as  material  constants,  independent  of  geometry, 
crack  size,  and  other  boundary  and  loading  conditions.    The  failure  conditions 
for  Mode  I  and  II  can  be  written  as: 


and 


F   .    a>/c~       >       K 


IC 


F   •    (  r  +       fa)     \AT     >      IC 


IC 


(2.49) 
(2.50) 


where 


a         =    applied  tensile  stress 
t  =    shear  stress 
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f        =   friction  coefficient  for  the  crack  surface 
F       =    a  geometrical  constant 

As  indicated  previously,  for  multi-modal  fracture,  the  Sih's 
strain  energy  density  factor,  Sq  ,  can  be  related  to  fracture  toughness  as: 

Sc    =    S    (KIC,  Knc,  Kmc)  (2.51) 

The  Sp  is  analogous  to  the  K^;  where  it  measures  the  resistance  of  material 
to  fracture  under  multi-modal  conditions.    Under  normal  tensile  stress 
conditions  where  Kq  =  Kjjj  =  0  , 

Sc    =    (1  -  2  y  )  (l  +  V   )        .      K     2  (2.52) 

2    77    E 

As  it  was  shown  previously,  similar  equations  can  be  also  developed  for 
other  fracture  modes. 
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B.  Fracture  Toughness  Measurement 

As  indicated  previously  the  plane  strain  fracture    toughness 
Kj£  is  a  material  constant  which  is  a  measure  of  the  material's  resistance 
to  fracture.    This  material  constant  has  been  shown  to  be  independent  of 
geometry    boundary  and  the  mode  of  load  application.    The  experimental 
procedures  for  fracture  toughness  measurement  and  evaluation  are  as 
follows: 

1.    Compliance  Calibration 

The  presense  of  flaws  and  cracks  within  a  body  affect  the  load-deflection 

characteristics  and  can  be  evaluated  using  the  material's  modulus  or  its 

inverse,  known  as  compliance.      Irwin  (L9)  has  shown  that  the  change  in  the 

spring  constant,  modulus  or  the  compliance  (inverse  of  modulus)  can  be  used 

as  a  measure  of  crack  extension  and  the  release  of  excess  elastic  energy. 

For  materials  tested  under  a  load  control  mode,  i.e. ,  "soft  systems,"  the 

energy  release  rate  is  written  as 

G    =  J_   -P2    •      3L  (2.53) 

2B  dc 

For,  stroke  control  testing,  i.e.,    stiff  machines  ,  the  energy  release  rate, 

however,  is  given  by: 

where  P,  A,  L,  B  are  load,  displacement,  compliance  and  specimen  width, 
respectively.    The  compliance  procedures  have  been  almost  exclusively 


57 


applied  to  Mode  I  problems   .    The  advantage  of  this  method  is  that  it  only 

requires  one  size  of  prototype  specimen.     Furthermore,  the  material  to  be 

tested  and  the  protype  need  not  be  the  same.    However,  it  is  required  that 

material  to  be  homogenous,  isotropic,  having  a  well  defined  modulus  of 

elasticity,    as  well  as  having  dimensions  proportional  to  the  real  specimens. 

For  analysis  purposes,  the  compliance  measurement  is  carried  out  on 

a    specimen  with  an  initial  notch  c0,  and  then  the  notch  is  extended  by  a 

small  increment  Ac  and  the  compliance  is  measured  again.    This  procedure 

is  continued  until  a  compliance-crack  length  relation  is  resulted.    The  rate  of 

change  of  compliance  with  crack  length  is  then  calculated  from  the  shape  of 

compliance-crack  length  as  given  by   5L    .    The  strain  energy  release  rate  G 

dc 

and  the  strain  intensity  factor  K  (K2  =  EG)  is  then  calculated  using  the  equa- 
tions given  previously.    The  comparison  of  the  Kj£  results  obtained  by 
experimental  compliance  method  and  analytical  procedures  have  been  shown 
to  be  in  excellent  agreement. 

2.    Analytical  Method 

In  this  method,  the  strain  analysis  procedures  are  used  to  derive 
mathematical  formulations  between  the  strain-intensity  or  strain-energy 
release  rate  and  the  crack  length  and  specimen  geometry.    One  such  proce- 
dure is  the  equation  given  by  Winne  and  Wundt,  as  shown  by 


°IC  =    (1  -  p2)    h   an     f(c/d) 
E 


(2.55) 
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As  it  was  shown  previously,  for  small  crack     f(c/d)   =  — -: —  (1  -  c/d)  , 
whereas  for  deeply-notched  beams,    f(c/d)   =    0.52. 

Other  analytical  equations  used  for  stress-intensity  factor  calculations 
were  given  by  equations  (2.37),  (2.40^  and  (2.41)  in  Chapter  I,  Part  II.    In 
addition,  numerous  approximate  analysis  procedures  have  also  been  recom- 
mended by  various  researchers  (2)  (Figure  2.14). 

Similarly,  the  use  of  the  polynomial  form  of  the  Ktq  equation  has  also  been 
reported  in  studies  of  asphaltic  concrete  (16)  and  Portland  cement  concrete  (20). 
One  such  equation  (2855)  used  for  Portland  cement  concrete  is  written  as: 

_     6M     f     -2c     ,.,.  /J4     "I  1/2 


[     TF     h(c/d)    J 


IC  Bd2    |_       n  J  <2-56> 

where  h(c/d)     =     10. 08(c/d)2    -    1.225 (c/d)    +    0.1917 

The  detailed  discussions  of  Ktq  measurements  for  pertinent  civil  engineering 
materials  will  be  discussed  in  a  separate  section.    In  all  such  analyses,  the 
maximum  load  Pq  shall  be  calculated  in  accordance  with  the  following  speci- 
fied procedures: 

(i)  Draw  an  initial  tangent  to  the  load-displacement  curve 

(line  OA)  as  shown  in  Figure  (2. 15). 

(ii)  Draw  the  second  line  O  P5  through  the  origin  with  slope 
less  than  the  slope  of  initial  tangent  OA.  P5  is  the  load 
at  intersection  of  O  P5  with  the  load-displacement 

59 


-!|N 

U 


CO 


CO 

•r-l 

H3 

S3 

.o. 

Cl> 

a 

+ 

00 

as 

o 

. 

a 

**■*• 

CO 

tj 

CM 

hn 

^^ 

a 

hf) 

.d. 

1    t-. 

-o 

XI 

C 

faD 

LO 

<M      'O 

o 

xi 

4— > 

0) 

a 

-u 

a 

•r-i 

T3 

■  1— ( 

t- 

T3     O 

CM 

PQ 

CD 

4-J 

CO 
■  1— t 

CD 

-Q 

■r-l 
O 

a 

1 

d 

0) 
O 

CM* 

1 

CO 
CJ 

r-l 

13.66  (c/ 
+    25.22  ( 

bX) 
C 
•r-* 

C 
CD 

00 

II 

II 

4_> 
C 

^ 

i 

o 

•i-1 

•  I-" 

CI) 

CO 

Sh 

CD 

O 

o 

rH 

r-l 

II 

a 

rH 

a 

a 

CD 

o 

u 

a 

3 

i 

i 

^H 

o 

cd 

Cl 

CO 

CO 

*5 

^ 

CO 

ft 

>H 

II 

r-t 

II 
CXI 

11 

CO 

II 

- 

ft  — 

1 

t- 

1 

CM 

"ft 

co       rr 

r*H          CM         > 

CNJ 

ft 

CO 

© 

LO 

© 
© 

CO 


TJ 


■* 


CM 


CO 


"tf 


CM 


CO 


CO 


CO 


CO 


CM 


CM  CM 

A 


CM 


CM 


CO 

a 

CD 

a 

•r-l 

o 
0 
a 
co 

S3 

PQ 

r-l 
O 

<*-! 

CQ 

o 

•1-4 

CvJ 
r- 

rQ 

u 


« 


■* 


CM 

CD 

r-l 

S) 

•I— I 

ft 


60 


3 


bo 

o 

-J-> 

s 

0) 

o 

i — i 

a 

CO 


03 

u 
o 
o 

0) 

« 

a 

CD 

s 

o 

I— I 

a 

CO 

•rH 

Q 

i 

o3 


o 

03 

CD 


03 

a 
o 

•1-H 

P. 


io 


0) 

u 


pBOl 


61 


curve.     To  determine  Pq,  if  the  load  at  every  point  on  the 
curve  which  precedes  P5  is  lower  than  P5,  then  Pq  =  P5 
(Type  I).    If,  however,  there  is  a  maximum  load  preceding 
P5,  which  exceeds  it,  then  this  maximum  load  is  equal  to 
Pq.     To  check  the  validity  of  Pq,   select  a  point  on  the  load- 
displacement  such  that  P  =  0.  8  P^.     Measure  the  distance  X 
along  the  horizontal  line  from  tangent  OA  to  the  curve.     Simi- 
larly, measure  the  horizontal  distance  Y  from  OA  to  P5.    If 
the  ratio  of  X  to  Y  is  greater  than  0.  25,  then  the  Kjc  is  not  a 
valid  test  and  strain- intensity  calculated  is  considered  as  a 
candid  value  of  strain- intensity  rather  than  as  a  true  measure- 
ment of  this  material  constant. 
3.  J-Integral  Criterion 

In  the  previous  sections,  the  failure  criterion  for  linear  elastic  mechan- 
ics with  small  scale  yielding  were  presented.    It  was  pointed  out  that  the 
strength  of  crack  tip  singularity  can  be  represented  by  the  stress  intensity  fac- 
tor, K,  and  its  lower  limiting  critical  value,  Kj^. 

For  large  scale  yielding,  Rice  (68,  32)  has  proposed  the  path- independent 
J- Integral  as  a  criterion  for  failure.     For  linear  elastic  material,  Jj£  is  identical 
with  Gjc  parameter  and  can  therefore  be  related  to  the  fracture  toughness,  Kj^. 
For  a  general  case,  however,  the  J- Integral  is  interpreted  as  potential  energy 
difference  between  two  identically  loaded  and  neighboring  cracks;    i.e. , 
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J    =    -  -? —  /  unit  thickness  (2.57) 

where  U    =    potential  energy 

c   =    crack  length. 

To  evaluate  the  fracture  criterion,  the  load-displacement  is  plotted  for  various 

crack  sizes  as  used  for  compliance  procedures.    At  a  given  total  deflection,  the 

area  under  the  curve  is  measured  using  the  planimeter.    The  energy  represented 

by  the  area  under  the  curves  are  plotted  versus  the  crack  length  as  shown  in 

Figure  2.16.    The  slopes  of  the  lines  are  given  by  Equation  (2.57).    It  should 

be  noted  that  load-displacement  procedures  used  should  measure  the  total  energy 

input  accurately.    The  use  of  crack  opening  displacement  method  for  J-integral 

analysis  is  not  satisfactory.     Furthermore,  the  deflection  at  failure  needs  to  be 

measured.    In  Figure  (2.16),  variations  of  U  vs.  c  per  unit  thickness  are  shown. 

The  slope  of  this  diagram,  i.e.  dU/dc  versus  deflection  is  shown  in  Figure  (2.17) 

The  experimental  observations  would  also  provide  an  estimate  of  deflection  at 

failure.    As  an  example,  considering  that    Af,  the  deflection  at  failure,  is  0.025 

inch  (Figure  2.18),  then  J  at  failure  can  easily  be  estimated. 

C.    Experimental  Results 

The  critical  toughness,  K„,  is  affected  by  the  specimen's  geometry  as 
well  as  the  boundary  and  loading  conditions.  In  order  to  obtain  a  reproducible 
value  for  the  lower  limiting  critical  toughness  KT„,  which  can  be  considered 
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as  a  material  constant,  the  state  of  stresses  should  closely  approximate  a  plane 
strain  condition.    This  requirement  and  restriction  on  the  stress  state  implies 
that  the  size  of  plastic  zone  must  be  small  as  in  comparison  with  the  dimensions 
of  the  specimen.    The  radius  of  the  plastic  zone  is  given  by  an  approximate 
relation: 

r    =    0.05   I— -=^1  (2.58) 


where  IC.-  is  fracture  toughness  and  <r    is  the  material  yield  strength  in  tension. 
To  meet  the  plane  strain  condition,  ASTM  has  established  a  minimum  require- 
ment for  the  plate  thickness  B  and  crack  length  c  as  given  by: 

/KIC\ 
B    >  2.5    I 1  (2.59a) 

and  /k- \  2 

and  specimen  depth  d  >2B.    (2.59b) 


■ » m 


The  deviations  from  such  a  size  requirement  will  result  in  higher  values  of 
Kjc  or  plane  stress  representation  of  fracture  toughness,  as  shown  in  Figure 
(2.19).    It  should  be  noted  that  as  the  specimen's  thickness  increases,  the 
stresses  acting  on  the  crack  front  also  increase,  leading  to  more  triaxial  state 
of  stresses  at  the  crack  front.    Such  a  condition  corresponds  to  a  plane  strain 
state.    In  Figure  (2.20),  a  schematic  view  of  the  effect  of  specimen  size  on 
fracture  toughness  and  the  test  results  on  metals  is  shown. 

Another  important  variable  is  the  mode  of  load  application  and  the  so- 
called  soft-hard  characteristics  of  the  testing  machine.    The  load-controlled 
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testing  systems  are  generally  considered  as  soft  machines.      A  typical  load- 
deformation  characteristic  of  specimens  tested  under  a  soft  machine  is  shown 
in  Figure  (2.21).    It  should  be  noted  on  the  machine-specimen  instability  path 
as  shown  in  Figure  (2.21),  a  plane  stress  failure  condition  might  occur.      On 
the  other  hand,  stroke-controlled  testing  systems  are  ideally  stiff  machines  in 
which  fracture  occurs  in  plane  strain  condition  (Figure  2.22). 

The  fracture  toughness  Kjc  has  also  been  shown  to  be  affected  by  tem- 
perature as  well  as  rate  of  load  application.      In  Figure  (2.2  3),  the  variation 
of  Kjq  with  temperature  is  shown.    It  is  also  noted  that  yield  strength     <ry 
is  also  influenced  by  temperature.      Generally,  at  high  temperatures,      cr 
decreases,  whereas  Ktq  exhibits  an  increase  with  temperature. 

The  experimental  data  on  metals  also  show  that  there  is  an  inter- 
relation between  tensile  strength  and  fracture  toughness,  Kjq.      It  is  gene- 
rally observed  that  Kjq   and    c-y   are  inversely  interrelated.      Although  there 
is  no  basic  reason  why  one  should  expect  a  simple  relation  between  Kjq  and  <r 
however,  such  a  relation  can  be  a  useful  approximation  for  Kjq  calculation. 

Experimental  data  on  metals  have  shown  such  a  relation  can  often  be  approxi- 

-3 
mated  by     Ktq  =    <r       .      Similarly,  an  empirical  formula  relating  Kjq  and 

o-y  has  been  presented  as  follows  (Figure  2. 24): 


'IC  3 

where  °"y    =     yield  strength 


KTn      =      |  t  •    <ry  •   n2-    *fl2  (2.60) 


n     =     strain  hardening  exponent   a-  =    C    £ 


a  constant 
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E    =     modulus  of  elasticity 
£f     =     true  strain  at  fracture. 

It  should  be  noted  that  all  material  would  not  necessarily  satisfy  such  an 
interrelation  between  a    and  Kjq. 

In  Tables  2. 1  and  2.  2,  the  available  Kjq  and  Gjq  data  on  numerous  engi- 
neering materials  are  shown.  The  fracture  toughness  data  on  concrete,  asphalt 
and  other  paving  materials  are  also  available  on  a  limited  basis. 

The  fracture  analysis  on  concrete  was  carried  out  in  the  early  sixties 
by  M.   F.  Kaplan  (37)  who  performed  tests  on  concrete  beams  with  crack  simu- 
lating notches.     The  critical  strain- energy-release  rate,  Gp,  associated  with 
the  rapid  extension  of  the  crack,  was  determined  by  two  methods:    the  analytic 
method  and  the  direct  experimental  method.      Similar  beams  with  different 
depths  of  notch  gave  Gc  values  which  are  in  close  agreement.      However,  Gq, 
values  of  beams  of  larger  dimensions  were  higher  than  those  of  smaller  dimen- 
sions.     Some  differences  in  Gq  values  were  reported;  the  discrepancy  is  possibly 
due  to  the  method  of  G^  determination.    Kaplan's  conclusion  was  that  the  critical 
strain-energy-release  rate  may  be  ascertained  by  suitable  analytical  and  experi- 
mental methods,  and  the  fracture  strength  of  concrete  containing  cracks  can 
thereby  be  predicted. 

Moavenzadeh  and  Kuguel  (45)  also  used  notched-beam  specimens  of 
cement  paste,  mortar  and  concrete  to  study  their  fracture  properties.     They 
found  out  that  the  fracture  work  of  the  paste  increases  by  the  introduction  of 
solid  particles.      This  is  due  to  the  multiplicity  of  the  crack  growth  during  the 
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TABLE  2.1   FRACTURE- TOUGHNESS  VALUES* 


Material 


Ultimate  Strength 


u 


ksi 


Critical  Stress- 
Intensity  Factor 
Kq-  ksi  Vin 


A517F  Steel  (AM) 


120 


170 


AISI  4130  Steel  (AM) 

AISI  4340  Steel  (VAR) 
AISI  4340  Steel  (VAR) 
AISI  4340  Steel  (VAR) 
AISI  4340  Steel  (VAR) 
AISI  4340  Steel  (VAR) 


170 

300 
280 
260 
240 
220 


100 

40 
40 
45 
60 
75 


300M  Steel  (VAR) 
300M  Steel  (VAR) 
300M  Steel  (VAR) 
300M  Steel  (VAR) 
300M  Steel  (VAR) 


300 
280 
260 
240 
220 


40 
40 
45 
60 
75 


D6AC  Steel  (VAR) 


240 


40-90 


H-ll  Steel  (VAR) 
H-ll  Steel  (VAR) 
H-ll  Steel  (VAR) 


320 
300 

280 


30 
40 
45 


12Ni-5Cr-3Mo  Steel  (VAR) 


190 


220 


18Ni  (300)  Maraging  Steel  (VAR) 
18Ni  (250)  Maraging  Steel  (VAR) 
18Ni  (200)  Maraging  Steel  (VAR) 
18Ni  (180)  Maraging  Steel  (VAR) 


290 
260 
210 
195 


50 

85 
120 
160 


9Ni-4Co-0.  3C  Steel  (VAR) 


260 


60 


Al  2014- 
Al  2024- 
Al  2219- 
Al  2618- 
Al  7001- 
Al  7075- 
Al  7079- 
Al  7178- 


T651 

T851 

T851 

T651 

T75 

T651 

T651 

T651 


70 
65 
66 
64 
90 
83 
78 
83 


23 
23 
33 
32 
25 
26 
29 
24 


*    These  values  are  taken  from  References  (73)(74)  and  (75)  and  are  to  be  considered 

as  nominal  values  only. 
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TABLE  2.2      CRITICAL  STRAIN  ENERGY  RELEASE  RATE  VALUES  * 

Material  Critical  Energy 

Release  Rate 
GIC -lb/in 

Dural  1.60  x  103 

Key  Steel  5.71  x  102 


Brass  3.43  x  102 


Teak  Wood  6.85  x  10 

Cast  Iron  4.57  x  10 

Cellulose  2.28x10 

Polystyrene  1. 14  x  10 

Polymethylmethacrylate  5.71 

Epoxide  Resin  3.77 

Polyester  Resin  2.51 

Graphite 

Alumina 

Magnesia 

Glass 


5.71  x  10"1    - 

•    1.14 

4.57  x  10"1 

1.14  x  10"1 

4.57  x  10"2 

*    Most  of  these  data  were  obtained  from  a  notched  bar  under  three-point 
bending   (73) 
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fracture  process  in  the  specimens.    The  extent  of  the  internal  cracking  was 
measured  using  a  quantitative  microscopy  technique.    The  result  of  such  crack 
measurements  indicated  that  the  true  fracture  work  of  concrete  (determined 
by  accounting  for  multiplicity  of  cracks  in  the  specimen)  is  lower  than  that 
of    cement  paste.    This  is  attributed  to  the  preference  of  cracks  to  propagate 
through  the  interface  of  paste  and  aggregate,  which  in  general  is  of  lower  bond 
strength  than  the  paste  matrix. 

Naus  and  Lott  (50)  defined  an  effective  fracture  toughness  for  concrete, 
assuming  that  it  is  a  homogeneous  material.    The  conclusions  of  his  studies 
were  as  follows:    The  effective  fracture  toughness  of  concrete  increases  with 
age  of  concrete,  increase  of  maximum  size  of  coarse  aggregate  and  increase 
of  gravel-cement  ratio.    There    is    no  apparent  effect  of  varying  the  water- 
cement  ratio  on  the  effective  fracture  toughness  of  the  concrete.    However,  a 
decrease  of  effective  fracture  toughness  was  observed  with  increase  of  air 
content  in  the  concrete. 

K.  P.  George  (10) has  also  studied  the  applicability  of  the  Griffith 
theory  of  brittle  fracture  to  the  fracture  of  soil-cement  material.    He  has 
investigated  the  possibility  of  using  the  concept  of  a  critical  strain- energy- 
release  rate  Gc,  or  the  fracture  toughness  Kc,  determined  from  Gc,  in 
predicting  the  crack  propagation  rate  in  soil-cement  base.    Critical  strain- 
energy-release  rate,  Gq,  and  the  fracture  toughness,  K^,  of  five  different 
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soil-cement  mixtures  were  experimentally  determined  on  notched  beams  and 
their  variations  with  respect  to  several  parameters  were  investigated.    Some 
relevant  conclusions  were  drawn: 

1.  The  value  of  the  critical  strain-energy-release  rate,  G^,  is  more 
or  less  independent  of  the  notch  depth. 

2.  Gq increases  as  the  clay  content  of  the  soil  becomes  greater. 

3.  Gq  increases  slightly  as  temperature  is  lowered  to  the  freezing 
point  of  water  and  thereafter  increases  rapidly  as  temperature 
drops  below  the  freezing  point. 

4.  Gq  increases  slightly  with  the  rate  of  loading. 

The  variations  of  Gq  according  to  soil  texture  and  the  rate  of  loading 
as  well  as  its  independence  of  notch  geometry  indicate  that  Gq  is  truly  a 
measure  of  physical  properties  of  the  materials.    He  also  examined  the  crack 
propagation  in  model  base  slabs  constructed  to  simulate  field  conditions.    It 
has  been  shown  that  the  critical  stress  intensity  factor,  Kq,  computed  from 
Gq,  tends  to  be  related  to  the  crack  propagation  rate. 

The  application  of  fracture  mechanics  principles  to  analysis  and  eva- 
luation of  asphaltic  systems  has  been  carried  out  by  Moavenzadeh  (44) , 

Herrin(l6),  Majidzadeh  et  al.   (33-41),    Monismith  (46),  and  Blight  (3). 
Moavenzadeh  investigated  the  application  of  Griffith  theory  to  fracture  analysis 
of  asphalt  cements  tested  at  low  temparatures.    Utilizing  the  Winne  and  Wundt 
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equation  presented  previously,  the  critical  strain- energy-release  rate  Gjq  was 
calculated  for  three  asphalt  cements.    E  was  concluded  that  at  sufficiently  low 
temperatures,  the  asphalt  cement  behaved  as  a  brittle,  amorphous  material 
satisfying  the  requirements  of  Griffith  theory.    It  was  pointed  out,  as  shown  in 
Figures  (2.25)  and  (2„26)  that  Gtq  for  asphalt  cement  is  a  rate  and  temperature 
dependent  parameter.      The  degree  of  aging  and  the  asphalt  type  also  influence 
the  critical  strain-energy-release  rate  (Figure  2.27). 

Majidzadeh,  Kauffmann  and  Ramsamooj  (35) (36) (37),  have  also  investi- 
gated the  fracture  characteristics  of  asphaltic  materials  at  various  tempera- 
tures and  rates  of  loading.      Kauffmann  (38)  studied  the  variations  of  fracture 
toughness  Kj£  with  rate  of  loading  and  temperature  (Figure  2.  28).      For  sand- 
asphalt  mixtures  tested  at  14,  73  and  32°F,  it  was  reported  that  the  critical 
stress  intensity  factor  increased  with  the  rate  of  load  application.      Similarly, 
the  fracture  toughness,  Kjq,  increased  with  the  test  temperatures.     A  series 
of  1"  x  1"  x  12"  beam  specimens  were  prenotched  with  notch  sizes  varying  from 
3/32"  to  12/32"  in  depth  before  testing.    The  experimental  results  indicated 
that  for  such  geometrical  conditions,  the  K™-.  is  independent  of  c/d  ratio  and 
remains  as  a  material  constant. 

Monism ith  and  Salam  (46),  in  their  investigation  of  asphalt  mixture 
behavior,  conducted  two  types  of  fracture  tests: 

(i)  Single- edge-notched  bend  specimens  of  1.5"  x  2.0" 

x  15"  tested  in  four-point  bending, 
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(ii)    Single-edge-notched  tension  tests  on  specimens  of  1.5  x  1.5  x 
4.5". 
To  maintain  a  constant  loading  condition,  all  specimens  were  tested  at  a  load- 
ing time  of  about  0.5  seconds.    The  experimental  results  presented  by 
Monismith  and  Salam(46)  indicate  that  at  low  temperatures  (less  than  10°F) 
KTC  increases  with  asphalt  content   up  to  8%  level.    At  higher  temperatures, 
the  measured  KTC  indicated  a  decrease  with  asphalt  content.    It  has  been 
postulated  that  an  increase  in  the  plastic  flow  and  crack  tip  blunting  might 
have  been  responsible  for  the  reduction  in  Kw-,  parameter. 

The  effect  of  aggregate  gradation,  asphalt  hardness  (Figs.  2.29-2.32), 
and  void  content  has  been  reported  to  be  of  lesser  importance  when  specimens 
are  tested  at  temperatures  of  40°F  or  higher. 

The  aggregate  type,  on  the  other  hand,  were  reported  to  have  a  signi- 
ficant influence  on  the  fracture  toughness. 

Blight  (  3  )  has  also  investigated  the  fracture  properties  of  paving 
materials.    Utilizing  the  modified  Griffith  equation,  the  effective  fracture 
surface- energy  y     is  given  by 

<r    =   constant  [E 7\2  (2.61) 


W 


where  y   is  the  effective  fracture  surface- energy  or  the  sum  of  the  true 

surface-energy    y,  and  the  plastic   work  ,        yp  .    The  experimental  results 
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presented  indicate  that  the  effective  fracture  surface-energy     y     of  asphaltic 
material  is  related  to  strain  at  failure  and  its  value  differs  for  various  mate- 
rials.   However,  it  is  noted  that  fracture  surface-energy  is  independent  of 
the  method  of  testing  employed. 
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HI.     MECHANISTIC  MODELS   FOR  FATIGUE 

A.  Review 

The  "fatigue"  is  considered  as  a  process  of  progressive  damage  or 
deterioration  under  repeated  or  cyclic  loading  which  eventually  leads  to 
failure  of  the  structural  system.    As  a  phenomenon  of  a  highly  complex 
nature,  it  involves  a  localized  progressive  structural  change  within  the 
material,  which  can  be  subdivided  into  three  stages:    crack  initiation,  crack 
growth  and  terminal  state  of  fracture.    The  occurrance  of  these  processes 
in  a  material  system  result  in  a  gradual  weakening  of  the  structural  compo- 
nent.   During  the  stage  of  crack  initiation,  microcracks  are  originated  at 
centers  of  inpurities,  flaws,  and  microstructural  defects.    These  centers  of 
strain-incompatibility,  when  subjected  to  reversed  cyclic  strain,  are  believed 
to  be  responsible  for  crack  initiation  process.    The  second  stage  of  fatigue 
process  is  the  crack  propagation  which  at  first  is  preceded  by  a  zone  during 
which  micro-macro  crack  transitions  occur. 

The  process  of  crack  propagation  and  the  terminal  state  of  fracture 
have  been  treated  by  numerous  theories.    The  introduction  of  fracture  mech- 
anics principle  into  analysis  of  fatigue  of  material  systems  has  provided  an 
analytical  method  of  classifying  the  crack  severity.    It  has  also  presented  a 
rational  scheme  for  the  life  expectancy  calculations  of  structural  systems. 

In  such  a  mechanistic  approach  it  is  postulated  that  the  crack  growth 
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is  a  consequence  of  the  changing  of  the  crack  tip  profile.    During  a  cyclic 
deformation,  a  crack  will  undergo  a  phenomenon  of  blunting  and  resharpening. 
In  tensile  loading  cycles,  the  crack  tip  tends  to  open  first  and  then  becomes 
blunted  as  the  plastic  zone  forms  and  spreads  ahead  of  the  crack  tip.    During 
the  unloading  cycle,  the  elastic  contraction  of  the  material  surrounding  the 
crack  imposes  a  residual  compressive  stress  on  the  plastically-deformed 
material  at  the  crack  tip.    This  reduces  ductility  and  resharpens  the  crack 
which  aids  the  growth  of  the  crack  in  the  next  loading  cycle.    This  process 
leads  to  slow  crack  growth  until  the  crack  reaches  a  critical  size,  where 
unstable  fracture  occurs.    The  factors  affecting  the  rate  of  crack  growth  are 
stress  or  strain  amplitude  and  the  defect  size  which  determine  the  stress 
intensity  at  the  crack  tip. 

The  mechanism  of  the  fatigue  crack  growth  can  be  explained  in  classi- 
cal terms  of  the  energy  balance  at  the  crack  tip.    The  work  of  the  external 
force  at  the  crack  tip  is  divided  into  stored  elastic  energy,  surface  energy 
required  to  form  cracks,  and  deformation  energy  required  for  some  irrever- 
sible structural  distortions.    The  rate  at  which  a  crack  may  propagate  and 
the  path  it  follows  depends  entirely  on  this  energy  balance  at  the  crack  tip. 
When  there  is  a  large  amount  of  plastic  deformation,  the  crack  may  become 
blunted  and  does  not  propagate  or  propagates  very  slowly.    On  the  other  hand, 
when  the  elastic  energy  released  exceeds  the  energy  demanded  for  creating 
new  surfaces,  the  crack  will  propagate  along  the  path  which  requires  the 
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minimum  amount  of  energy  to  create  new  surfaces. 

The  mechanistic  formulation  of  the  fatigue  process  is  credited  to  Paris, 
Gomez  and  Anderson  (59)  who  first  introduced  the  application  of  the  stress- 
intensity  factor  to  fatigue  crack  propagation  rates.    In  1963,  Paris  and  Erdogan 

dc 
found  from  experimental  data  that  the  crack  propagation  rate,  —  ,  was  propor- 
tional to  the  fourth  power  of  ^  K  for  a  number  of  materials.    This  law  of 
crack  growth  is  expressed  as  —  =  A  Kn,  where  A  and  n  are  material  con- 
stant (n  =  4.0). 

The  fourth  power  relation  has  been  justified  by  consideration  of  the 
energy  absorption  within  the  entire  plastic  zone  ahead  of  the  crack  tip  ( 6  ) . 

However,  experimental  data  available  on  many  engineering  materials  indicate 

dc 
that  n  might  be  smaller  than  four.    In  Figures  (2.  34  -  2.  36),  typical     —  vs.K 

relations  for  metals  are  provided. As  it  is  noted,  the  power  constant  n  varies 
over  a  wide  range  depending  on  material  characteristics  and  testing  condi- 
tions. 

The  data  available  for  asphaltic  mixtures  also  indicate  that  constant 
n  varies  depending  on  asphaltic  mixture  characteristics.    For  fine  grained 
asphaltic  mixtures,  as  sand  asphalt  and  surface  course  asphaltic  concrete 
with  small  size  aggregate,  a  more  rapid  crack  propagation  has  been  reported 
(33),  (34),  (35) .    The  n  =  4  condition  has  been  also  reported  by  Majidzadeh,  et.  al. , 
(40),  (36),  (37),  as  valid  for  sand-asphalt  mixtures  tested  under  monotonic  load- 
ing conditions.    The  random  or  block  loading  condition  also  has  been  found  to 
be  affecting  the  crack  propagation  process. 
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To  incorporate  the  effect  of  other  test  variables,  such  as  frequency, 
mean  load  and  load  range,  a  number  of  empirical  equations  have  also  been 
proposed  in  the  literature.     As  an  example,  Mukherejee  (49)  has  presented 
the  following  equation: 

dN     ~  '  <Kmax)  (*»K)  (2.63) 

where  f  is  frequency,  Kmax  and    £K  are  maximum  stress  intensity  factor  and 
Kmax  -  Kmin  respectively.     This  equation  indicates  that  frequency  of  loading 
influences  the  crack  growth  rate.     The  effect  of  rest  period  has  also  been  in- 
vestigated by  a  number  of  researchers.    It  has  been  argued  that  the  mechanism 
of  crack  propagation  under  dynamic  loading  and  static  conditions  are  similar 
except  for  the  rest  period  effect. 

The  experimental  results  by  Kauffmann  (28)  indicate  that  the  fatigue 
crack  growth  process  in  asphaltic  mixtures  is  considered  as  a  stress-independent, 
stress- interaction  type  phenomena.    It  has  been  observed  that  for  the  range  of 
variables  investigated,  the  rest  period  does  not  significantly  affect  the  fatigue 
life.     The  available  data  also  indicate  that  load  sequence  has  a  significant  effect 
on  the  fatigue  life.     The  delay  in  the  rate  of  crack  propagation  due  to  loads  of 
variable  amplitude  in  the  loading  sequence  is  attributed  to  the  beneficial  effects 
of  residual  stresses  at  the  crack  tip.    Although  the  mathematical  formulation 
of  the  delay  phenomena  has  not  yet  been  fully  developed,  the  reduction  in  the 
crack  growth  rate  can  be  evaluated  experimentally  (Figures  2.  37  and  2.38 ). 

The  crack  growth  model  has  been,  in  recent  years,  subject  to  numerous 
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modifications.    One  such  analysis  is  the  work  of  Cherepanov  (6)  who  utilized  the 

plastic  work  at  the  crack  tip  to  develop  a  relationship  between  total  plastic 

work  and  the  stress-intensity  factor.    This  development  is  presented  by  an 

equation  given  by:  2 

/'k2M  -  K2min>    ♦    In/  KC    -  K2max    V\ 

=  /*      —tt  I  rjt :  K2     j J     <2-64> 

min        / 


d  N  r        \  K(f  \    Kc"  -  Kt 


where     &   is  a  material  constant  given  by: 

P   =    a4        £_£c  (2#65) 

2    (T  z 

y 

In  this  equation,  a.  is  a  dimensionless  parameter  related  to  <r„/E  and  1/ 

where  ©■„  is  the  yield  stress  of  the  material. 

The  crack  growth  models  are  in  general  incorporated  into  three  load- 
ing variables:    mean  stress-intensity  factor  Kmax,  cyclic  stress-intensity 
factor  AK  =  Kmax  -  Kmin,  and  stress  ratio  given  by  R  =  K^^/K^^. 

The  effect  of  mean  load  and  load  range  has  also  been  incorporated  in 

Pearson  (62)  and  Porter    (63)  models  of  crack  propagation.    The  Pearson  (62) 

model  is  expressed  in  terms  of  R,  load  ratio  and  critical  stress  intensity 

factor.    This  model  can  be  written  as: 

dc    =  A(AK)n 

dN  (1  -  R)  Kc-  A  K  (2*65) 

K 

where  R  =        min    ,    and  K^is  the  critical  stress  intensity  factor.     Pearson 

Kmax 
(62)  has  indicated  that  for  materials  with  relatively  high  critical  stress  inten- 
sity factor,  the  load  ratio  R  does  not  significantly  affect  the  crack  growth 
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rate.    However,  for  conditions  where  Kmax  approaches  KjC ,  the  crack  growth 
data  should  be  modified  as  follows: 


dc       _       AKc(jK)n 

dN      ~      a-H)Kc-/lKm  (2'66) 


where  m  and  n  are  constants. 

In  a  more  recent  study,  Majidzadeh  et  al.   (41)  have  used  the  results  of 
linear  elastic  fracture  mechanics  to  study  the  fatigue  of  field  asphaltic  speci- 
mens.    The  analysis  of  fatigue  and  fracture  of  field  specimens  has  indicated  that 
the  power  n  may  not  always  be  4  and  might  vary  with  material  characteristics  and 
test  conditions.    There  are  indications  that  the  conditions  of  the  testing  and  the 
properties  of  the  asphaltic  mixtures  may  influence  this  power  law.    As  an  example, 

in  this  study,  a  power  law  of  the  second  degree  in  K  has  proved  to  be  of  the  high- 
est correlation. 

This  observed"  difference  in  the  exponent  of  the  power  law  suggests  that  it 
might  be  desirable  to  review  the  models  and  recent  advancements  in  the  theory 
underlying  the  crack  propagation  process.    This  was  done  in  the  hope  that  the 
variables  such  as  the  power  and  material  constants  could  be  related  to  the  funda- 
mental properties  of  asphaltic  concrete  and  the  testing  procedure. 

It  is  well  known  that  asphaltic  concrete  displays  a  delayed  response  due 
to  loading.      This  rheological  response  depends  on  the  environment  and  the 
rate  of  loading.    It  could  be  linear  or  non-linear,  depending  upon  those  condi- 
tions and  the  intensity  of  the  load.    Because  of  the  complexity  of  this  problem, 
linearity  of  the  material  will  be  assumed.    It  is  realized  that  this  assumption 
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may  not  hold  under  actual  conditions,  but  with  the  present  state  of  knowledge 
it  is  necessary. 

This  literature  review  will  not  be  concerned  with  the  micro-n±,:"hanics 
aspect  of  fracture  or  the  fatigue  problem.    It  will  concentrate  on  the  conti- 
nuum mechanics  perspective.    Asphaltic  concrete  is  assumed  to  be  linearly 
viscoelastic,  so  that  concepts  from  linear  elastic  fracture  mechanics  can  be 
used. 

It  is  not  possible  to  discuss  in  detail  the  vast  volume  of  research  that  has 
accumulated  in  recent  years  concerning  the  crack  propagation  problem. 
The  work  of  some  researchers  will  be  reviewed  briefly  in  preparation  for 
certain  applications. 
B.  Viscoelastic  Fatigue  Damage  Modeling 

The  review  of  literature  indicates  that  theoretical  approach  to  analysis 
of  generalized  crack  propagation  problems  involve  consideration  to  three 
stages.    The  first  stage  is  to  obtain  a  solution  for  the  boundary  value  problem, 
using  the  material  constitutive  equation  which  is,  in  this  case,  considered  as 
a  linearly  viscoelastic  material.    The  second  stage  is  the  use  of  a  failure  cri- 
terion and  ,  lastly,  a  crack  model  assumption  is  needed. 

1.    Viscoelastic  Boundary- Value  Problem 

Usually  the  solutions  to  viscoelastic  boundary- value  problems  can  be 
obtained  by  applying  the  classical  correspondence  principle.    This  is  con- 
venient as  the  solutions  to  elastic  boundary  value  problems  are  readily 
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available  and  can  be  transformed  to  obtain  the  solution  to  the  viscoelastic 
ease.    However,  it  should  be  noted  that  this  simple  correspondence  principle 
cannot  be  applied  to  fracture  problems  since  these  problems  involve  moving 
boundaries  (cracks).    Graham  (11)  established  a  useful  extension  for  crack 
problems.    This  extended  correspondence  principle  applies  with  certain 
restrictions,  as:    the  crack  size  should  increase  monotonically  with  time; 
the  stresses  are  independent  of  elastic  constants  and  finally  these  elastic 
constants  appear  in  a  separate  factor  in  the  displacement  expression.    The 
crack  closure  case  does  not  satisfy  the  first  restriction.      All  three 
restrictions  are  satisfied  for  all  quasi-static  problems  if  tractions  are  spe- 
cified on  the  boundaries,  tractions  on  crack  surfaces  are  self- equilibrating 
and  there  are  no  body  forces. 

To  obtain  the  solution  to  the  viscoelastic  BVP,  the  Laplace  transform 
of  the  elastic  solution  for  the  same  BVP  is  taken.    Next,  the  elastic  constants 
are  replaced  by  the  transformed  viscoelastic  expressions.     Finally,  the 
expressions  are  inverted  back  to  the  time  domain.    In  the  case  of  a  crack  with 
the  above  restrictions,  the  stresses  will  be  the  same  for  the  elastic  and  the 
viscoelastic  cases.    The  displacement,  however,  will  be  rather  complicated 
for  the  viscoelastic  case  and  will  involve  convolution  integrals.    The  com- 
plexity of  the  displacement  expression  is  the  reason  for  the  lengthy  mathe- 
matics and  some  of  the  simplifying  assumptions  involved. 

2.     Fracture  Criteria 

The  fracture  criteria  used  by  most  investigators  are  modifications  of 
the  Griffith  energy  balance  formulation.    Consider  the  energy  balance  equation: 
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U    =    V    +    T    +    D 


(2.67) 


where  U  =  the  work  done  by  the  external  loads 
V  =  the  elastic  stored  energy 
T  =  the  kinetic  energy 

D  =  irreversible  energy  such  as  surface  energy,  plastic  work, 
and  viscous  dissipation. 

The  dot  indicates  differentiation  with  respect  to  time.    Theoretically,  these 

energies  can  be  evaluated  for  the  particular  boundary  value  problem.    They 

can  be  used  in  conjunction  with  the  equation  of  motion  to  obtain  a  solution  for 

the  rate  of  crack  extension  (Erdogan,  8_  ).    In  practice,  this  is  very  difficult 

to  apply.    Simplifications  in  analysis  lead  to  the  assumption  of  a  plane 

symmetric  problem  with  no  body  forces,  and  neglecting  kinetic  energy.    The 

energy  balance  will  then  become: 


J7V,ds  =fffanW*      +   47T6 


(2.68) 


*f  being  the  specific  fracture  energy. 

To  simplify  this  quasi-static  problem  considerably,  a  local  energy 
balance  can  be  used  in  a  region  surrounding  the  crack  tip.    The  assumption 
implied  here  is  that  the  size  of  the  plastic  zone  is  small  compared  to  the 
region  where  the  energy  balance  is  applied.    Consider  the  following  two 
figures  (Knauss,   29  ): 

,A 


R 
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x  =  0  —J 


For  Figure  1,  the  energy  balance  is: 

U    =    V   +    D   +    2    S  (2.69) 

D  =  rate  of  change  of  dissipated  energy;  S  =  rate  of  change  of  surface  energy. 
For  Figure  2,  there  is  no  S  since  there  is  no  crack: 

4U    =    |V    +       y*Tn,  Uy    dx   +    ID  (2.70) 

Comparing  the  two,  it  is  evident  that: 

^        Uy     dx   =    S   =     Tf      a  (2.71) 

This  is  for  a  traction-free  crack.    This  local  energy  principle  applies  to  many 
BVP's  as  was  used  by  various  investigators  (29, 31,47 , 69,7_9) .    It  can  be  stated 
simply  as  with  crack  advancement  the  work  of  the  released  tractions  goes  into 
creating  new  surfaces. 

3.    Crack  Model  Assumption 

The  crack  models  presented  in  the  literature  are  primarily  those  of 
Griffith,  Barenblatt,  Dugdale,  etc.    The  Griffith  model  is  the  oldest  and 
simplest.      It  considers  the  crack  to  be  in  a  perfectly  brittle  material.     There 
is  no  yielding  at  the  crack  tip  although  the  stresses  can  be  very  high.    The 
resulting  equation  is: 


i 
o-  a2     =     K  (2.72) 
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where    o-     is  the  stress  applied  at  the  boundaries  of  the  infinite  plate,  a  is 
half  the  crack  length  and  K  is  the  stress  intensity  factor. 

The  Dugdale  model  assumes  a  thin  plastic  zone  at  the  tip  of  the  crack. 
This  special  configuration  yields  an  estimate  of  the  plastic  zone  size  as: 


<L      =    2  sin2      Jl_  _£_  *-2     o-2 

4     Y  8Y2 


(for      °"     <<1)  (2.73) 


where  at   is  the  size  of  the  plastic  zone,  Y  is  the  yield  stress,      <r  and  a  as 
above. 

The  Barenblatt  model  (4)  proposes  the  addition  of  cohesive  forces 
along  the  faces  of  the  crack  that  resist  the  opening  of  the  crack.    This  model 
is  advantageous  in  getting  rid  of  the  singular  solution  at  the  tip  of  the  crack 
inherent  in  the  above  models: 

NQ    =V~q"    ^m1!  (2.74) 

7T 

where  N     =    _I_  ,    L  is  some  integral  depending  upon  the  assumed  stress  dis- 
tribution,      <r      is  the  maximum  stress  in  the  cohesive  zone.    The  Barenblatt 
model  includes  the  Dugdale  model  as  a  special  case. 

It  is  possible  now  to  combine  the  results  of  the  boundary- value  problem 
solution,  the  crack  model,  and  the  failure  criterion  to  obtain  a  solution  for  the 
crack  problem  in  viscoelastic  media.    This  can  be  quite  involved  due  to  the 
presence  of  convolution  integrals  in  the  expressions  written  for  the  displace- 
ment.   To  get  rid  of  the  complexities  of  this  problem,  Williams  (76)  studied 
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the  simple  geometry  of  a  spherical  cavity  in  a  body  under  all-around  tension. 
Since  variations  in  the  crack  shape  result  only  in  the  change  of  a  constant 
in  the  criticality  condition,  the  study  of  a  simple  geometry  will  still  reflect 
the  viscoelastic  character  of  fracture.    In  this  manner  Williams  did  not  have 
to  assume  a  crack  model,  he  used  a  simple  cavity  for  which  the  BVP  solution 
is  simple.    In  this  manner  the  calculations  become  simpler.    He  applied  a 
global  energy  balance,  and  for  the  case  of  a  step  load,  he  obtained: 

^oc  =   4     /        T/a0 


1  -    (ao)  3  3  >      2    J(t0)  -  Jff 

b  S 


V- 


(2.75) 


where  aQ   =    cavity  radius 

b      =    the  radius  of  body 
T     =    the  fracture  energy 
J     =    the  creep  compliance 
Jg   =    the  glassy  compliance 

The  Griffith  solution  can  be  retrieved  from  the  above  formula  for 
the  glassy  case  by  putting  J  (tQ)  =  J 

Other  investigators  assumed  actual  crack  models  rather  than  a  simple 
cavity.    The  work  of  these  investigators  will  be  grouped  according  to  the 
models  they  assumed  and  are  presented  in  the  following  sections  of  this 
chapter. 
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4.    Analysis  Procedure 

a.    Griffith  Geometry 

Knauss  (30)  studied  the  Griffith  geometry  for  the  case  of  a  linear 
viscoelastic  material.    Using  a  local  energy  balance  equation  and  assuming 
the  stress  distribution  at  the  top  of  the  crack  to  be  given  by  the  linearly 
viscoelastic  solution,  he  obtained: 


V 


/  1^L  (2.76) 


oc  a 

J  (T) 

a 

where  a  (t)    =    half  of  crack  length 

a    =    the  size  of  the  crack  failure  zone. 
It  is  noted  that  this  equation  is  more  complicated  than  Williams'  result  because 
of  the  more  complex  geometry.    Since  the  crack  length  depends  on  time  and  its 
first  derivative  is  in  the  argument  of  the  creep  function,  the  resulting  equation 
is  a  nonlinear  first  order  differential  equation.    Its  solution  depends  on  some 
approximations  and  whatever  creep  function  is  assumed. 

Mueller  (47)  studied  the  problem  of  a  large  crack  in  an  infinitely  long 
strip  under  a  prescribed  strain  in  the  lateral  direction.    Superposition  of 
two  simple  problems  was  used  to  obtain  a  nondimensional  stress  intensity 
factor  and  the  size  of  the  crack  opening.    The  solution  was  obtained  assuming 
a  stationary  crack  with  a  sudden  application  of  the  prescribed  strain.    This 
same  problem  was  extended  to  study  crack  propagation  by  applying  a  local 
energy  balance  (48).    It  was  found  that  the  crack  tip  velocity  depends  on  geo- 
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metry,  applied  load  and  its  history,  and  material  properties  as  follows; 


i 


*oEr   =  J*      T/b  (2-77) 

G(f) 

a 


where  (      =    the  prescribed  step  strain 

Er    =    the  long  time  relaxation  modulus 
b     =    half  width  of  the  strip 
G     =    some  function  of  the  creep  compliance 
If     f0Er  is  thought  of  as  the  stress,  the  similarity  of  this  result  to  the  previous 
ones  is  immediately  obvious. 

For  crack  propagation  under  variable  load  histories,  Knauss  and 
Dietmann  (31)  used  the  local  energy  balance  approach,  assuming  the  dis- 
tribution of  the  stresses  ahead  of  the  crack  to  be  of  a  singular  form  and  of 
a  certain  variation.    The  displacement  corresponding  to  this  stress  distri- 
bution is  obtained  "by  using  the  extended  correspondence  principle.    The  power 
equation  is  applied  and  a  differential  equation  for  the  crack  tip  velocity  as  a 
function  of  the  time  history  of  the  stress  intensity  factor  is  obtained.    The 
equation  is  a  complicated  nonlinear  differential  equation  and  for  specific 
stress  histories,  might  be  simplified.    If  the  stress  intensity  factor  does  not 
change  rapidly,  then  the  result  is  the  same  as  in  Reference  (48).    The  case  for 

a  cyclic  strain  input  was  also  examined  and  some  graphic  results  were  given. 
It  is  pointed  out  that  the  complexity  of  the  differential  equation  makes  it  of  little 

practical  interest  in  its  present  form. 
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It  is  noted  that  in  all  of  the  above  work,  no  special  consideration  of 
the  character  of  the  material  at  the  crack  tip  was  made.    This  zone  is  usually 
very  highly  stressed;  it  is  probably  in  a  plastic  state  or  might  even  be  in  a 
disintegrated  state  though  not  yet  cracked.    Although  the  behavior  of  the 
material  in  this  zone  is  not  well  understood,  one  might  be  able  to  make 
some  assumptions  that  will  simplify  the  calculations.    The  material  around 
the  crack  tip  is  probably  highly  nonlinear  and  therefore  making  such  assump- 
tions will  not  necessarily  lead  to  less  exact  results  than  the  linearity  assump- 
tion that  is  implied  in  the  previously  discussed  literature. 

It  has  been  shown  by  some  investigators  that  the  material  at  the  tip  of 
the  crack  is  crazed  and  of  a  lower  density  than  the  rest  of  the  material.    The 
material  behavior  here  is  very  similar  to  metal  yield  and  is  in  a  wedge-like 
shape.    Wnuk  and  Knauss  (81)  used  these  observations  to  assume  the  visco- 
plastic  model  of  Crochet  for  the  tip  of  the  crack  while  the  bulk  material  is 
linearly  viscoelastic.    A  solution  by  Graham  (11)  for  a  crack  opened  by  a 
normal  pressure  acting  on  its  surface  was  used  along  with  the  assumed  yield 
model.    With  some  assumptions,  the  final  nonlinear  differential  equation  will 
be: 


y2    -    y'y    =    y2        / P  (t)  y    +    Q  (t)\  (2.78) 


where  y(t)  =  Y(t)  -  A,  Y(t)  is  the  time  dependent  yield  and  A  is  a  material  con- 
stant from  Crochet  model. 

\/~o    x  c 
P(t)    = - i//(t)  where    X  =  1/2  (1  -  v  )  (1  +     v),  C  is  a  material 

g 
constant  from  the  Crochet  model,  E„  is  glassy  modulus,  and  0  (t)  is  the  slope 
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of  nondimensional  creep  compliance  J(t)/J(o)  and  finally  Q(t)  =  2^ ^-^ —  (A+B)  kjj 

\  2         Eg 

where  B  is  also  a  material  constant  in  the  Crochet  model.    This  differential 
equation  is  also  very  difficult  to  solve  in  general.    A  simple  case  of  a  Maxwell 
model  is  solved  in  (81).    Another  possibility  is  to  assume  a  time  independent 
yield  and  the  result  will  be  simpler  (81) . 

b.    Dugdale  Crack  Model 

This  model  was  used  by  many  workers  in  the  field  of  fracture  because 
of  its  simplicity.    Wnuk  (78),  (79),  (80)  used  it  to  simplify  his  mathematical 
derivations.    Wnuk  proposed  a  nonlinear  differential  equation  to  describe  the 
quasi- static  subcritical  growth  of  a  crack.    He  derived  this  equation  in  a 
simple  way  by  applying  a  local  energy  balance: 

(a)  Amount  of  energy  supplied  during  a    6  L  advance  of  the  crack  is: 

(2.79) 

The  first  term  is  for  creation  of  new  surfaces  and  the  second  is  for  strain 
redistribution. 

(b)  Amount  of  energy  demand  is :  <pc     8L 

6  U    \        +       6U    )  dp  =0  /2.80) 

.'.       6L  6a  /         dL  ,      c 

O  L 

This  equation  can  be  put  in  a  simpler  form: 

G   +    M    =    Gc  (2.81) 

where  G   =  the  energy  release  rate 

M  ■'=  the  slow  growth  operator 

Gq  =    the  critical  energy  for  onset  of  rapid  fracture. 
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AU=      MM     SL   +    6_U  JO         5    L 


For  the  case  of  an  elastic  solid  with  a  plastic  region  surrounding  the 

crack  tip  the  G  and  M  operators  are  (78) : 

6  /*a 

G   =    2  Y     ufcl      +    2    Y     ^j-  /  u    (x)    dx  and 


do-  6 

M  =  2   Y      df       67 


/'  • 


then 


M   =   2    Y      clo       6  (  A    u     v 

dl     ^o^  V    3       W 

Utip   =    ~§f  and      A      =        g^-        G 


(2.82) 


where  Y  =  the  constant  yield  stress  in  the  plastic  zone 

or  =  the  applied  stress 

1  =  half  the  crack  ie.igth 

a  =  half  the  crack  length  plus  the  yield  zone  length 

u(x)  =  the  displacement  normal  to  the  crack  plane. 

For  the  viscoelastic  case  with  a  plastic  yield  zone  the  differential 
equation  can  be  extended  as  simply  (78): 

(G   +    M)   vj/  t      j-     J     =    G  (2.83^ 

where  xjf   =  the  normalized  creep  compliance  function,  and 

A    -  the  inherent  opening  distances. 

This  differential  equation  can  be  simplified  by  (78) : 

(a)    Assume  linear  range        w«l  and,  therefore,        —  <<  1, 

■  L 

G  =  2  Y  V 


(2.84) 
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and  so  M  becomes: 


M     = 


7TEG 

12  n  Y5 


dG 
d  o- 


do- 
dl 


where    H     = 


l-i/- 


plane  stress 


plane  strain 


(2.85) 


and  the  differential  equation  becomes: 


1    + 


7T21(T  do- 


6Y' 


dl 


G       \J/ 


A 


1 


=     GC 


(2.86) 


(b)  Assume  a  Dugdale  model  for  the  crack  tip  shape.    Using  a  solution 

for  the  displacement  normal  to  the  crack,  Wnuk  finally  puts  the  differential  equa- 


tion for  the  elastic  case  in  the  form: 
d_£    _  3    2-C& 


2      CW 


for  plane  crack 


dX    _  3   2-CX2 
dC    '  2      £2A3 
where      Q  _    *  o  \  '      ° 

2y      '  y 

• 

For  the  viscoelastic  case,  the  differential  equation  is: 

.2 


for  penny- shaped  crack 


(2.87) 


(!<* 


d£ 


2-tf 


'Vx  +  c\dx       2-Cx2 


plane  crack 
penny-shaped  crack 


(2.88) 


where 


C  =  ¥ 


t=0 


I 


Though  these  differential  equations  have  no  closed  form  solutions,  it  will 
be  seen  later  how  they  can  be  of  help  in  determining  the  form  of  the  crack  propa- 
gation law  with  some  approximations. 
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c.  Barenblatt  Crack  Model 

In  the  above  Dugdale  model  and  the  previous  modified  Griffith  model, 
assumptions  were  made  regarding  the  failure  zone  behavior.    In  the  Barenblatt 
equilibrium  model,  the  behavior  of  the  failure  zone  is  arbitrary.     This  is  much 
more  general  and  covers  the  others  as  special  cases. 

Schapery  (69)  utilized  the  Barenblatt  equilibrium  model,  a  local  energy 
criterion,  and  the  elastic  solutions  of  Williams  and  Barenblatt.     By  superposing 
the  elastic  solution  of  Williams  for  a  crack  in  a  plane  strain  due  to  external  load- 
ing and  the  solution  of  Barenblatt  for  a  crack  with  a  stress  of   cr    in  the  failure 

f 

zone  with  no  external  loading,  one  obtains  an  expression  for  the  stress  intensity 
factor  and  the  displacement  (69): 

N0     =      V*l     «TM    l±  (2.89) 


where  L     =        /        f(qr) )  dr) 


1     =         q  and 

°"f 
f       =    — —     is  the  distribution  of   <rf  in  terms  of  maximum  stress. 

*m 
N0  is  the  stress  intensity  factor  of  Barenblatt  and  is: 

KI 
NQ     =     "^ZT  (2.90) 


V27T 


The  expression  for  N    can  be  thought  of  as  a  relation  for  the  size  and  nature  of  the 
failure  zone.     The  displacement  expression  is: 
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2Ce  ^m  ^'^ 

J^T       —qT        I2         £  H(f)       ,  where  (2.91> 


— = — ^ ^    ,    plane  strain 

Ce     =        J  ft  - 

4/E         ,    plane  stress 


i  f1  _1I     dn 

2       o*^     dri   Vn~ 


and         H     =     Heaviside  function. 

Now  that  the  elastic  stress  and  displacement  expressions  are  available,  the 
extended  correspondence  principle  can  be  utilized  to  get  the  viscoelastic  counter- 
parts.    Here,  an  assumption  of  a  constant  Poisson's  ratio  will  simplify  matters 
and  only  the  uniaxial  creep  will  enter  the  equations. 

When  the  displacement  expression  is  transformed  to  the  viscoelastic  case, 
a  convolution  integral  arises.     The  assumption  of  small  curvature  on  a  log-log  plot 
of  the  creep  compliance  is  used  to  reduce  the  convolution  integral  for  displacement 
to  a  simple  product  form  by  using  a  power  law  representation  for  the  creep  com- 
pliance.    The  final  step  is  applying  the  local  energy  balance  and  Schapery  obtains: 

Cv  fa)    =       7rN2  where  (2.92) 

o 

~  ~  J-/n       a 

C     (ta)  is  the  creep  compliance  with  the  argument,    t„       =       A  — —    t  and  is 

some  function  of   n   which  arises  in  using  the  power  creep  curve;   a    is  the  failure 

zone  size  and  can  be  found  from: 

a-       *?No 

2     2 

°"m  I 
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r    is  the  fracture  energy,  NQ  is  the  stress  intensity  factor  of  Barenblatt,  and  for 

KI 

the  opening  mode  is  equal  to       ~/=^r~    >      °"      is  the  maximum  of  the  failure  zone 

stress  distribution  and  L  is  some  integral  which  depends  on  the  assumed  stress  dis- 
tribution,    r,  (T     and  L  are  considered  fracture  properties  and  experiments  should 
help  in  their  determination. 

Knauss  (29)  also  used  the  Barenblatt  equilibrium  model  because  of  its 
generality  and  the  finiteness  of  stresses  at  the  crack  tip.    He  points  out  that  the 
Dugdale  Model  is  a  special  case  of  the  Barenblatt  model.    He  uses  two  criteria 
for  crack  propagation:    the  energy  criterion  and  the  maximum  strain  criterion. 
His  final  results  have  a  convolution  integral  in  the  expressions  written  for  the 
crack  propagation  rate.     It  is  pointed  out  that  this  work  is  very  similar  to  the 
above  work  by  Schapery,  except  that  no  approximation  was  used  by  Knauss  to 
eliminate  the  convolution  integral. 

5.  Theoretical  Models  For  Fatigue  Process 

The  next  logical  step  is  to  extend  the  above  concepts  and  methods  of  solution 
to  problems  involving  repeated  loads.     This  is  not  as  straightforward  as  it  might 
appear  and  further  mathematical  simplification  and  assumptions  will  be  needed. 
Only  the  procedures  which  result  in  an  explicit  relation  for  the  crack  propagation 
law  will  be  considered  here.     The  others  which  can  be  solved  by  other  techniques 
but  for  which  an  explicit  form  is  not  readily  available  will  not  be  considered  further. 

Wnuk  (80)  uses  his  resulting  differential  equation  to  study  the  case  for  re- 
peated loads.     Fatigue  is  viewed  as  a  sequence  of  extensions.    Consider  the  case 
for  the  plane  crack: 
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,2     Tft     A      C    ,    dP  2-^P2 

(i     *P     +     "P    >    dT        =  j^  (2-93) 

Assuming  crack  length  constant  in  a  cycle,  it  can  be  shown  that: 

_d£_     =     J_        ,2  /    4  4      \  c^  3  3 

dn  12  \  Pmax  P0      J  6<P>         (P  P0    > 

(2.94) 
If    (3    is  assumed  to  be  above  the  minimum  cycle  stress  and  <  ($>  =   2f        , 

where  f  is  the  frequency,  the  final  result  will  be  after  changing  to  dimensional 

quantities: 


t-M^y^ 


The  first  term  being  the  familiar  Paris  law,  the  second  term  is  due  to  the 

rate  dependency  of  the  material.    Another  similar  result  was  obtained  by  Wnuk 

but  with  slightly  different  assumptions,  such  as  the  possibility  of  extension  during 

unloading:  /  a7  \A  „2  U  K  \  J  A  K  \*  1 

dl/dn=-^-    4a^TJ+9Cf(irJ|  (2.96) 

where     a=j£*  ^=kV*{£} 

Schapery's  (69)  analytical  approach  can  be  solved  to  get  an  explicit  relation 


for  the  tip  velocity: 


4  r 


or 


cv  <*«  )      =     ~^W  (2-97> 

v  o 


c   «    x1/n   4L,   =   *H2 

v         n  a  ttn2 


Putting  in  the  power  creep  law,     C  (t)     =     C^  t    : 
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vl/n       a      n  4f 

ci(  An     -T-)    =  irss 


0 


an  4f  (2'98) 


Ci/n    ~7*      =  ^Ng 


2 


but 


Hi 


a2  if 


2    2  \n 
"o    J        77-Np 


*  Nn  \        ttNo 


ci  An  I  ^Tf/'  4r      "  a 

x     m    1  ' 

/  2n+l  X1/11  2+  2/n 

so  da     _       /  C-,  An    «"  )  N0 

dt  -  (— - —  ;   -*v-  (2.99) 

m    1 


For  crack  growth  in  one  cycle,  it  is  assumed  that  crack  growth  per  cycle 

is  small  and  it  is  furthermore  assumed  a  wave  shape  of  stress  intensity  factor, 

No 
w    =    .77—      .      N        is  the  maximum  value  of  the  stress  intensity  factor.     Then, 
Norn  Om 

da  „     2(l4) 

d¥    =     BN0m  <2-100> 


where  /  2n+l         \  1/n  t+L, 


V'n~   2  T2  dt 


m    1 


Schapery  argues  that  if  bitumen  completely  surrounds  the  tip  and  for  tempera- 
tures not  close  to  the  glass  transition  temperature,  the  slope  may  be  taken  as  1  and 
the  Paris  power  law  is  retrieved.    Usually  though,     0  4    n    4.5  for  viscoelastic 
materials,  implying  the  power  is    6. 
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C.  Generalized  Fatigue  Models:    Dimensional  Analysis 

The  above  review  of  fracture  and  fatigue  in  viscoelastic  media  indicates  that 
there  is  much  similarity  in  the  fatigue  process  for  elastic  and  viscoelastic  media. 
The  next  few  pages  will  deal  with  fatigue  laws  obtained  from  dimensional  analysis, 
applying  in  general  to  any  material. 

Kang  and  Liu  (26),  assuming  a  plane  strain  condition  and  small  scale  yield- 
ing, applied  dimensional  analysis  to  the  problem  of  crack  growth.     For  these  assump- 
tions, he  argued,  the  crack  tip  region  can  be  scaled  by  the  size  of  the  plastic  zone. 
This  means  that  the  stresses  and  strains  at  the  geometrically  similar  points  are  the 
same  for  various  levels  of  stress  intensity  factors.    Therefore,  stresses  and  strains 
must  be  the  same  within  crack  increments    da.  if  d  a  is  proportional  to  the  size  of 
the  plastic  zone.    If  the  material  is  homogeneous,  then   da/dN  is  proportional  to  r  • 

da/dN     =      fx(R)  rp  ,  (2.101) 

2 
where    r    is  the  size  of  the  plastic  region  and  is  proportional  to  (AK)    .    Therefore, 

da/dN      =      f2(R)  (      K)2  (2.102) 

Li  analyzing  actual  data,  Kang  and  Liu  found  that  the  power  of  the  above  model 
varied  depending  on  the  level  of   Ak.     The  following  figure  illustrates  this  variation. 
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Figure  2.39    Typical  fatigue  crack  propagation  data  for 

2024-T351  aluminum  alloy  (from  Kang  and  Liu) 
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The  figure  shows  that  for  low    A  K,  the  power  (the  slope)  is  4  or  more. 
Also,  for  high    AK,  the  slope  is  also  4  or  more.     For  the  intermediate  region, 
the  slope  is  close  to  2.    Kang  and  Liu  also  point  out  that  these  regions  might  be 
better  fitted  by  curves  rather  than  straight  lines.    This  has  been  the  recent  ex- 
perience at  Ohio  State  University  as  well  (41). 

To  illustrate  this  point  further,  the  dimensional  analysis  of  Cherepanov 
and  Halmanov  will  be  reviewed.  Assuming  continuous  crack  growth  and  a  con- 
stant dissipation  energy,  they  obtain: 


9  2  2  2 

—    =      -      BfKmax   "    Kmin     +   in     KC     ~     Kmax  (2  103\ 

*C  KC    ~     Kmin, 


where    |$  is  a  material  constant  that  depends  on  yield  stress,     <r     ,  modulus  of 
elasticity,   E,  Poisson's  ratio     v    ,  prehistory  of  loading       <r0,  and  the  fracture 
toughness    K^  .     In  general, 

P  =     a^,-^,    *    )    lA  <2-104> 

2  <rz 


2  2 


being  a  non-dimensional  constant. 

/Kc  -    KTTiaxl 
To  simplify  the  above  rate  of  propagation  model,  the  second  term   ln|-^ j, 


\K"C  -     K 


mm/ 


can  be  expanded  in  a  Taylor  series  around  (0,  0) .    Substituting  back  in  the  above 
model  yields: 


da 

dN 


ft  I  ^ax   "   Kmin  Kmax   "    Kmin  \ 

(2.105) 
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This  is  not  quite  similar  to  Liu's  result  but  it  confirms  the  presence  of  high 
order  powers.    If   K^     =    0  ,  then  one  obtains: 

6  \ 


da      _     a  I    KmoY  K 


max         ,         nmax 


+ 


dN  V      Kc  Kc 


+    .   .   .      J  (2.106) 


It  is  obvious  here  that  for  different  ranges  of  Kmax,  one  obtains  different 
powers  if  only  one  term  is  used  to  describe  the  process.     This  was  the  case  in  the 
work  done  by  Majidzadeh  et  al.   (41)  where  the  generality  of  the  power  law  was  in- 
vestigated. 

If  Kmax/K£  is  high  (low  cycle  fatigue),  then  higher  order  terms  should  be 
retained  and  the  last  part  of  Figure  2. 39 is  obtained.      Here  the  points  may  not  lie 
on  a  straight  line  as  was  observed.      If  Kmax/K^is  low  (high  cycle  fatigue),  then 
higher  order  terms  can  be  neglected  as  compared  to  the  leading  term  and  the  Paris 
power  law  is  obtained. 

It  is  conceivable  to  think  of  a  region  between  the  above  two  extremes  where 
a  second  power  law  might  replace  the  above  model.    In  this  region,  the  above  law 
will  still  describe  the  process  better.     But  for  the  sake  of  simplicity,  if  one  wants 
to  choose  only  one  power,  then  the  power  would  be  two.    In  a  previous  work  by  the 
authors,  this  was  found  to  be  precisely  the  case.    Asphaltic  concrete  specimens  were 
fatigued  at  different  temperatures  but  approximately  in  the  same  range  of  crack 
growth  rate  (41).    It  was  found  that  the  best  fit  would  be  obtained  if  one  uses  three 
or  more  terms  and  the  next  best  would  be  a  two  term  fit.    The  one  term  fit  led  in- 
variably to  a  power  of  2.     This  one  term  fit  was  considered  good  because  of  the 
simple  power  law  that  results. 


123 


part  n 

LIST  OF  REFERENCES 


1.  ASTM  Special  Committee  on  Fracture  Testing  of  High-Strength  Metallic 

Materials,  "Progress  in  Measuring  Fracture  Toughness  and  Using 
Fracture  Mechanics,"  Materials  Research  and  Standards,  Vol.  4, 
No.   3,  pp.  107-119,  March  1964. 

2.  ASTM  Special  Technical  Publication  463,  Review  of  Developments  in 

Plane  Strain  Fracture  Toughness  Testing,  National  Aeronautics 
and  Space  Administration,  September  1970. 

3.  ASTM  Special  Technical  Publication  381,  Fracture  Toughness  Testing 

and  Its  Applications,  1965. 

4.  Barenblatt,  G.  I. ,  "The  Mathematical  Theory  of  Equilibrium  Cracks 

in  Brittle  Fracture,"  Advances  in  Applied  Mechanics,  Vol.  VII, 
Academic  Press,  1962. 

5.  Brown,  W.  F. ,  Jr.  and  Srawley,  J.  E.,  Plane  Strain  Crack  Toughness 

Testing  of  High-Strength  Matallic  Materials,  ASTM  STP  410, 
American  Society  for  Testing  and  Materials,  1967. 

6.  Cherepanov,"~G.  P.  and  Halmanov,  H. ,  "On  The  Theory  of  Fatigue 

Crack  Growth,"  Engineering  Fracture  Mechanics,  Vol.  4,  1972. 

7.  Deverall,  L.  L.  and  Lindsey,  G.  H. ,  "A  Comparison  of  Numerical 

Methods  for  Determining  Stress-Intensity  Factors,"  ICRPG  Paper, 
1969. 

8.  Erdogan,  F. ,  "Crack  Propagation  Theories,  "  Fracture,  An  Advanced 

Treatise,  Vol.  2,  (H,  Liebowitz,  Editor),  Academic  Press,  New 
York,  New  York,  1971. 

9.  Fukakura,  Juichi,  "The  Effect  of  Specimen  Thickness  and  Temperature 

on  Plastic  Yielding  and  Fracture  Properties  of  Mild  Steel, "  Engineer- 
ing Fracture  Mechanics,  Vol.  6,  pp.  231-244,  Pergamon  Press,  1974. 

10.    George,  K.   P. ,  "Theory  of  Brittle  Fracture  Applied  to  Soil  Cement," 
Journal  of  the  Soil  Mechanics  and  Foundation  Division,  ASCE,  Vol. 
96,  No.  SM3,  pp.  991-1010,  May  1970. 


124 


11.  Graham,  G.  A.  C. ,  "The  Correspondence  Principle,  of  Linear  Viscoelasticity 

Theory  for  Mixed  Boundary  Value  Problems  Involving  Time  Dependent 
Regions,"  Quarterly  of  Applied  Mathematics,  26,  1968. 

12.  Griffith,  A.  A. ,  "The  Phenomena  of  Rupture  and  Flow  in  Solids,"  Philosophical 

Transactions  of  the  Royal  Society  of  London,  221,  163-168,  1921 

13.         ,  "The  Theory  of  Rupture,"  Proceedings  of  the  First  Inter- 
national Congress  in  Applied  Mechanics,  1924. 

14.  Gross,  B.  and  Srawley,  J.  E. ,  "Stress- Intensity  Factors  for  Three  Point 

Bend  Specimens  by  Boundary  Collocation,"  NASA  TN  D-3092,  Lewis  Re- 
search Center,  Cleveland,  Ohio. 

15.  Hahn,  G.  T.  and  Rosenfield,  A.R.,  "Sources  of  Fracture  Toughness:    The 

Relation  Between  Kjc  and  Ordinary  Tensile  Properties  of  Metals, "  Appli- 
cations Related  Phenomena  in  Titanium  Alloys,  ASTM  STP  432,  American 
Society  for  Testing  and  Materials,  pp.  3-31,  March  1968. 

16.  Herrin,  M.  and  Bahgat,  A.G. ,  "Brittle  Fracture  of  Asphaltic  Mixes,"  Pre- 

sented at  the  Annual  Meeting  of  tne  Association  of  Asphalt  Paving  Tech- 
nologists,  February,  1968. 

17.  Hoeppner,  D.  W„  and  Knopp,  W.  E. ,  "Prediction  of  Component  Life  by  Appli- 

cation of  Fatigue  Crack  Growth  Knowledge, "  Engineering  Fracture 
Mechanics,  Vol.  6,  1974. 

18.  Irwin,  G.  R. ,  "Analysis  of  Stresses  and  Strains  Near  the  End  of  a  Crack  Tra- 

versing a  Plate,  "  Jj:_of_A^pl^Iechajiics,  Trans.  ASME,  79:361-4,  1957. 

19.  Irwin,  G.  R.  and  Kies,  J.  A.,  "Critical  Energy  Rate  Analysis  of  Fracture 

Strength,"  Welding  Journal  Research  Supplement,  Vol.  33,  pp.  193s  - 
198s,  1954. 

20.  Irwin,  G.  R. ,  "Fracture,"  Encyclopedia  of  Physics,  Vol.  6,  pp.  557-590, 

Springer  Press,  Berlin,  1958. 

22.  Irwin,  G.R.  and  Wells,  A.  A. ,  "A  Continuum  Mechanics  View  of  Crack  Pro- 

pagation, "  ^tajlu£gJcal_Review,  10  (58),  pp.  223-270,    1965. 

21.  Irwin,  G.R.  and  Kies,  J.,  "Fracturing  and  Fracture  Dynamics,"  Welding 

Journal  Research  Supplement,   February  1952. 

23.  Irwin,  G.  R. ,  Structural  Mechanics,  Pergamon  Press,  Great  Britain, 

p.  557,  1960. 


125 


24.  Irwin,  G.  R.,  "Fracture  Mechanics, "  Proceedings,  First  Symposium  on 

Naval  Structures  Mechanics,  MacMillan  Press,  1958. 

25.  Irwin,  G.  R.,  "Plastic  Zone  Near  a  Crack  Tip  and  Fracture  Toughness," 

Seventh  Sagamore  Ordinance  Materials  Research  Conference,  August 
1960.   . 

26.  Kang,  T.  S.  and  Liu,  H.  W. ,  "Fatigue  Crack  Propagation  and  Cyclic 

Deformation  at  a  Crack  Tip, "  International  Journal  of  Fracture 
Mechanics,  Vol.  10,  1974. 

27.  Kaplan,  M.  F. ,  "Crack  Propagation  and  the  Fracture  of  Concrete," 

Journal  of  the  American  Concrete  Institute,  Vol.  58,  No.  5,  pp.  591- 
610,  1961. 

28.  Kauffmann,  E,  M. ,  "The  Application  of  Fracture  Mechanics  Concepts 

to  Predict  Cracking  of  Flexible  Pavements,"  Ph.D.  Dissertation, 
The  Ohio  State  University,  1973. 

29.  Knauss,  W.  G. ,  "On  the  Steady  Propagation  of  a  Crack  in  a  Viscoelastic 

Sheet:  Experiments  and  Analysis, "  Deformation  and  Fracture  of 
High  Polymers,  H.  H.  Kausch,  J.  A.  Hassell,  and  R.  I.  Jaffee, 
Editors,  Plenum  Press,  1974. 

30.  Knauss,  W.  G. ,  "Delayed  Failure — The  Griffith  Problem  for  Linearly 

Viscoelastic  Materials,"  International  Journal  of  Fracture  Mechanics, 
Vol.  6,  1970. 

31.  Knauss,  W.  G.  and  Dietmann,  H. ,  "Crack  Propagation  Under  Variable 

Load  Histories  in  Linearly  Viscoelastic  Solids,"  International  Journal 
of  Engineering  Science,  Vol.  8,  1970. 

32.  Kobayashi,  Albert  S. ,  Editor,  Experimental  Techniques  in  Fracture 

Mechanics,  The  Iowa  State  University  Press,  Ames,  Iowa,  1973. 

33.  Majidzadeh,  K.;  Ramsamooj,  D.  V.  and  Chan,  A.  T.,  "Fatigue  and 

Fracture  of  a  Bituminous  Paving  Mixtures,"  AAPT,  1969. 

34.  Majidzadeh,  K.;  Kauffmann,  E.  M.  and  Ramsamooj,  D.  V.,  "Application 

of  Fracture  Mechanics  in  the  Analysis  of  Pavement  Fatigue,"  AAPT, 
Vol.  40,  1971. 


126 


35.  Majidzadeh,  K.  and  Ramsamooj,  D.  V.,  "Development  of  Testing  Pro- 

cedures and  a  Method  to  Predict  Fatigue  Failures  of  Asphaltic  Con- 
crete Pavement  Systems,"  Final  Report,  Project  RF  2873,  The  Ohio 
State  University  Research  Foundation,  March  1971. 

36.  Majidzadeh,  K. ;  Kauffmann,  E.  M.  and  Saraf,  C.  L. ,  "Analysis  of 

Fatigue  of  Paving  Mixtures  from  the  Fracture  Mechanics  Viewpoint," 
Fatigue  of  Compacted  Bituminous  Aggregate  Mixtures,  ASTM  STP 
508,  1972. 

37.  Majidzadeh,  K. ;  Ramsamooj,  D.  V.  and  Chan.  A.,  "Analysis  of  Fatigue 

and  Fracture  of  Bituminous  Paving  Mixtures— Phase  I,  RF  2845,  The 
Ohio  State  University  Research  Foundation,  May  1970. 

38.  Majidzadeh,  K.  and  Kauffmann,  E.,  "Analysis  of  Fatigue  and  Fracture 

of  Bituminous  Paving  Mixtures— Phase  II, "  The  Ohio  State  University 
Research  Foundation,  RF  2845,  March  1971. 

39.  Majidzadeh,  K.  and  Ramsamooj ,  D.  V.,  "Development  of  Testing  Pro- 

cedures and  Method  to  Predict  Fatigue  of  Asphalt  Concrete  Pavement 
Systems,"  RF  2873,  The  Ohio  State  Research  Foundation,  March  1971. 

40.  Majidzadeh,  K.  and  Kauffmann,  E.,  "Verification  of  Fracture  Mechanics 

Concepts  to  Predict  Cracking  of  Flexible  Pavements,"  RF3200,  The 
Ohio  State  University  Research  Foundation,  June  1973. 

41.  Majidzadeh,  K.;  lives,  George;  Makdisi-Ilyas,  F.  and  Saraf,  C.  L. , 

"A  Laboratory  and  Field  Durability  Study  of  Asphaltic  Mixtures," 
Engineering  Experiment  Station  (EES  356),  The  Ohio  State  University, 
1974. 

42.  Markstrom,  Karl,  "On  Fracture  Toughness  and  Its  Size  Dependence  for 

Steels  Showing  Thickness  Delamination, "  Engineering  Fracture 
Mechanics,  Vol.  4,  pp.  593-603,  Pergamon  Press,  Great  Britain,  1972. 

43.  McElivy,  A.  J.;  Boettner,  R.  C.  and  Johnson,  T.  L. ,  "On  Foundation  and 

Growth  of  Fracture  Cracks  in  Polymers,"  Fatigue,  An  Interdisciplinary 
Approach,  Syracuse  University  Press,  1964. 

44.  Moavenzadeh,  Fred,  "Asphalt  Fracture,"  Professional  Paper  P67-5, 

Massachusetts  Institute  of  Technology,  Cambridge,  Massachusetts, 
February  1967. 


127 


45.  Moavenzadeh,  F.  and  Kuguel,  R.,  "Fracture  of  Concrete,"  Journal  of 

Materials,  Vol.  41,  No.  3,  pp.  497-579,  1969. 

46.  Monismith,  C.  L.;  Hicks,  R.  G.  and  Salam,  Y.M.,  "Basic  Properties 

of  Pavement  Components,"  Report  FHWA-RD-72-19,  Federal  Highway 
Administration,  University  of  California  at  Berkeley,  September  1971. 

47.  Mueller,  H.  K. ,  "Stress  Intensity  Factor  and  Crack  Opening  for  a  Linearly 

Viscoelastic  Strip  with  a  Slowly  Propagating  Central  Crack,"  Inter- 
national Journal  of  Fracture  Mechanics,  Vol.  7,  No.  2,  1971. 

48.  Mueller,  H.  K.  and  Knauss,  W.  G. ,  "Crack  Propagation  in  a  Linearly 

Viscoelastic  Strip,"  Journal  of  Applied  Mechanics,  Transactions, 
ASME,  June  1971. 

49.  Mukherejee,  B.  and  Burns,  D.  J.,  "Growth  of  Part  Through  Thickness 

Fatigue  Cracks  in  Sheet  Polymethylmethacrylate,"  Engineering  Frac- 
ture Mechanics,  Vol.  4,  No.  4,  December  1972. 

50.  Naus,  D.  J.  and  Lott,  J.  L.,  "Fracture  Toughness  of  Portland  Cement 

Concretes,"  Proceedings,  ACI,  Vol.  66,  pp.  481-489,  1969. 

51.  Neuber,  H. ,  "Theory  of  North  Stresses,  "  Translation  published  by  J.  W. 

Edwards  Co.,  Ann  Arbor,  Michigan,  1946. 

52.  Orowan,  E. ,  "Energy  Criteria  of  Fracture,"  Welding  Journal  Research 

Supplement,  pp.  157-160,  March  1955. 

53.  Orowan,  E.  and  Felbeck,  D.  K. ,  "Experiments  on  Brittle  Fracture  of 

Steel  Plates,"  Welding  Journal  Research  Supplement,  pp.  570-575, 
November  1955. 

54.  Orowan,  E.,  Fatigue  and  Fracture  of  Metals,  John  Wiley  &  Sons,  Inc. 

1952. 

55.  Orowan,  E.,  "Fundamentals  of  Brittle  Behavior  of  Metals,"  Fatigue  and 

Fracture  of  Metals,  John  Wiley  &  Sons,  Inc.,  139-167,  1952. 

56.  Paris,  P.  C,  "The  Fracture  Mechanics  Approach  to  Fatigue,"  Pro- 

ceedings of  the  10th  Sagamore  Conference,  Syracuse  University  Press, 
107-132,  1965. 

57.  Paris,  P.C.  and  Si h,  G.  C.,  "Stress  Analysis  of  Cracks,"  Fracture 

Toughness  Testing  and  Its  Applications,  ASTM  STP  381,  39-76,  1965. 


128 


58.  Paris,  P.  C,  Fatigue— An  Interdisciplinary  Approach,  Syracuse  Uni- 

versity Press,  1964. 

59.  Paris,  P.  C.;  Gomez,  M.  P.  and  Anderson,  W.  E.,  "A  Rational  Analytic 

Theory  of  Fatigue,"  Trend  in  Engineering  (University  of  Washington), 
Seattle,  Washington,  Vol.  13,  No.  1,  1961. 

60.  Paris,  P.C. ,  "The  Growth  of  Cracks  Due  to  Variation  in  Load,"  Ph.D. 

Dissertation,  Lehigh  University,  September  1962. 

61.  Paris,  P.C.  and  Erdogan,  F.  J.,  "A  Critical  Analysis  of  Crack  Propa- 

gation Laws,"  Jmirj2idj}f_F^asjx^  Trans.  ASME,  Series  D, 

Vol.  85,  pp.  528-553,   1963. 

62.  Pearson,  S. ,  "The  Effect  of  Mean  Stress  on  Fatigue  Crack  Propagation 

in  Half  Inch  Thick  Specimens  of  Aluminum  Alloys  of  High  and  Low 
Fracture  Toughness,"  Engineering  Fracture  Mechanics,  Vol.  4, 
p.  9,  1972. 

63.  Porter,  T.  R.,  "Method  of  Analysis  and  Prediction  for  Variable 

Amplitude  Fatigue  Crack  Growth,"  Engineering  Fracture  Mechanics, 
Vol.  4,  No.  4,  1972. 

64.  Ramsamooj,  D.  V. ,  "The  Design  and  Analysis  of  the  Flexibility  of 

Pavements,"  Ph.D.  Dissertation,  The  Ohio  State  University,  1970. 

65.  Ramsamooj,  D.  V.;  Majidzadeh,  K.  and  Kauffmann,  E.  M. ,  "The 

Analysis  and  Design  of  the  Flexibility  of  Pavements,"  Proceedings, 
Third  International  Conference  on  the  Structural  Design  of  Asphalt 
Pavements,  University  of  Michigan,  1972. 

66.  Rice,  J.  R.  and  Rosengren,  G.  F. ,  Journal  of  the  Mechanics  and  Physics 

of  Solids,  Vol.  16,  pp.  1-12,  1968. 

67.  Rice,  J.  R. ,  Journal  of  Applied  Mechanics,  "Tranactions  of  the  American 

Society  of  Mechanical  Engineers,"  June  1968,  pp.  379-386. 

68.  Rice,  J.  R. ,  "Mathematical  Analysis  in  the  Mechanics  of  Fracture," 

Fracture,  H.  Liebowitz,  Editor,  pp.  191-311,  1968. 

69.  Schapery,  R.  A. ,  "A  Theory  of  Crack  Growth  in  Viscoelastic  Media," 

Mechanics  and  Materials  Research  Center,  Texas  A&M  University, 
College  Station,  Texas,  1973. 


129 


70.  Sih,  G.  C. ,  "A  Review  of  the  Three-Dimensional  Stress  Problem  for  a 

Cracked  Plate,"  International  Journal  of  Fracture  Mechanics,  Vol.  7, 
No.  1,  pp.  39,  1971. 

71.  Sih,  G.  C.  and  Liebowitz,  H„ ,  "Mathematical  Theories  of  Brittle  Fracture," 

Fracture,  H.  Liebowitz,  Editor,  pp.  67-190,  1968. 

72.  Sih,  G.  C,  "Three-Dimensional  Stress-State  in  a  Cracked  Plate, "  Pro- 

ceedings of  the  Air  Force  Conference  on  Fatigue  and  Fracture  of  Air- 
craft Structures  and  Materials,  September  1970. 

73.  Sih,  G.  C. ,  "A  Special  Theory  of  Crack  Propagation,"  Methods  of 

Analysis  and  Solutions  to  Crack  Problems,  edited  by  G.  C.  Sih, 
Wolters-Noordhoff  Publishing,  1972. 

74.  Srawley,  J.  E.;  Jones,  M.  H.  and  Brown,  W.  F. ,  "Determination  of 

Plane  Strain  Fracture  Toughness, "  Materials  Research  and  Standards, 
Vol.  7,  No.  6,  pp.  262-266,  June  1967." 

75.  Srawley,  J.  E. ,  "Plane  Strain  Fracture  Toughness,"  Fracture— An  Advanced 

Treatise,  edited  by  H.  Liebowitz,  Vol.  IV,  Chapter  II,  pp.  45-68,  1969. 

76.  Williams,  M.  L.,  "Initiation  and  Growth  of  Vise oelas tic  Fracture,  "  Inter- 

national Journal  of  Fracture  Mechanics,  1965. 

77.  Winne,  D.  HT  and  Wundt,  B.  M. ,  "Application  of  the  Griffith- Irwin  Theory 

of  Crack  Propagation  to  the  Bursting  Behavior  of  Disce,  Including 
Analytical  and  Experimental  Studies,"  Trans.  ASME,  Vol.  80,  pp.  1643- 
1655,  1968. 

78.  Wnuk,  M.  P.,  "Subcritical  Growth  of  Fracture  (Inelastic  Fatigue),  Inte r- 

national  Journal  of  Fracture  Mechanics,  1971. 

79.  Wnuk,  M.  P. ,  "Slow  Growth  of  Cracks  in  a  Rate  Sensitive  Tresca  Solid," 

Engineering  Fracture  Mechanics,  5,  1973. 

80.  Wnuk,  M.  P.,  "Prior  to  Failure  Extension  of  Flaws  Under  Monotonic 

and  Pulsating  Loadings,"  Engineering  Fracture  Mechanics,  5,  1973. 

81.  Wnuk,  M.  P.  and  Knauss,  W.  G. ,  "Delayed  Fracture  in  Viscoelastic 

Plastic  Solids,"  International  Journal  of  Solids  and  Structures,  6, 
1970. 


130 


PART  HI.      A  MATHEMATICAL  INVESTIGATION  OF  GEOMETRICAL 

MODELING  AND  FRACTURE  AND  FATIGUE  ANALYSIS  OF 
ASPHALTIC  BEAMS  ON  ELASTIC  FOUNDATION 


I.      INTRODUCTION 

Previous  research  studies  at  The  Ohio  State  University  (12,  L3,  L|, 
15),  have  led  to  the  formulation  of  a  fracture  mechanics  method  of  pavement 
fatigue  design  subsystem.      This  design  principle,   which  provides  a  unique 
methodology  for  fatigue  life  prediction  of  pavement  structures,  also  includes 
recommendations  for  material  testing  and  evaluation  procedures.      The  pro- 
posed laboratory  testing  method  for  fatigue  and  fracture  analysis  involves 
experimentation  on  beams  on  elastic  foundations.      In  such  a  fatigue  experi- 
ment, a  two-dimensional  elastic  system  of  asphaltic  beams  on  elastic  foun- 
dation is  subjected  to  repeated  loading,  and  it  is  generally  assumed  that 
the  fatigue  crack  propagation  rate  follows  a  power  law  relating  the  crack 

growth  rate  to  the  stress  -  intensity  factor,  K (dc/dN    =     AK  ).    The 

material  constants  associated  with  the  crack  growth  rate,  A  and  n,  are 
affected  by  such  parameters  as  load  frequency,  external  boundary  condi- 
tions, temperature,  and  statistical  distribution  of  flaws  in  the  materials, 
as  well  as  the  dimensions  of  the  models  used. 

In  such  an  analysis,  a  pavement  system  is  treated  as  a  two-dimen- 
sional elasticity  problem,  and  the  laboratory  investigation  of  beams  on  an 
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elastic  foundation  subjected  to  repeated  loading  is  carried  out  to  simulate  the 
dynamic  behavior  of  the  pavement  system  under  actual  traffic  loads.      How- 
ever, before  this  simulation  can  be  accurately  carried  out,  the  effect  of  the 
model  dimension  on  the  material  constants  needs  to  be  determined.    This 
report,  therefore,  is  concerned  with  the  development  of  proper  criteria  for 
the  design  of  such  models  and  will  include  procedures  to  establish  proper 
criteria  for  selection  of  suitable  beam  geometries,  determination  of  allow- 
able loads  for  fatigue  testing  experiments,  and  graphical  determination  of 
stress- intensity  calculations  for  beams  on  elastic  geometry  as  needed  for 
fatigue  and  fracture  analysis. 
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II.      A   REVIEW  OF  THEORETICAL   CONCEPTS: 


BEAM    ON  ELASTIC    FOUNDATION 


A.      An  Infinitely  Long  Beam  Resting  on  an  Elastic  Foundation  of  Finite  Depth 


To  investigate  the  behavior  of  an  asphaltic  beam  tested  on  an  elastic 
foundation,  first  consider  the  effect  of  a  sinusoidal  load  of  value  q  per  unit 
length  acting  directly  on  top  of  an  infinitely  long  elastic  foundation  of  unit 
width  and  a  finite  depth,  df,  resting  on  a  rigid  base  as  shown  in  Figure  1. 


q  =  qo 


cos  Ax 


I« 


Elastic  Foundation 


Ef,  "f 

'///A//////*/////// 


Rigid  Base 


Figure   3. 1 


Under  this  condition,  the  value  of  the  load  is: 


q(x)    =     qQ  cos   \  x 


(3.1) 


This  problem  is  a  two-dimensional  elasticity  problem,  the  governing  differen- 
tial equation  of  which,  in  terms  of  the  Airy's  stress  function  and  neglecting 
the  body  force,  is  given  by: 


v4* 


o     , 


(3.2) 
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where  <p  is  the  Airy's  stress  function  related  to  the  stress  components 
by: 


a 

X 

d2<P 
dy2 

°y    z 

2 

-   a  <p 

(3.3) 


T, 


a2  0 


xy       ax  ay 


and 


*74        a4      2    s4  a4 

Ox  ox  Oy  dy 


Besides  the  well-known  assumptions  of  linear  elasticity,  it  will  be 
assumed  that  there  is  no  shearing  stress  acting  at  the  interface  of  the  founda- 
tion  and  the  rigid  base.     With  this  additional  assumption,  the  boundary  con- 
ditions of  the  problem  are: 


T        =    0  \ 

XY  I 

>  aty    =    0  (3.5) 

o       -    -q  •  cos  Ax       ; 

y         H 


Txy    =   ° 

aty    =  df  (3.6) 

v      -    0 


where    v    is  the  vertical  displacement  along  the  y-axis.      It  can  be  shown 
that  the  stress  function  given  by: 
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r  2 

0(x,y)     =     cosXx     jc-^coshA.y    +    C£X  sinhXy    +   CgA(2sinhXy 

t— 

+    XycoshXy)    +   c4X(2coshXy   +   XysinhXy)  (3.7) 

satisfies  the  biharmonic  equation  (  3.  2).      Consequently,  the  corresponding 
stress  components  are: 


[2  2 

c-^  X  coshXy    +   c2  X  sinhXy    +   csX(2  sinhXy 

+    XycoshXy)    +   C4X(2coshXy    +   XysinhXy) 


X  cosXx 


c^y  sinhXy 


c-,coshXy    +   C2SinhXy    +   Coy  coshXy 


(3.8) 


(3.9) 


=     XsinXx     c-^AsinhXy    +   C2XcoshXy    +   Co(coshXy    + 

XysinhXy)    +   C4(sinhXy    +   XycoshXy)         .  (3.10) 

Applying  the  boundary  conditions  (  3„  5)  to  (  3.  9)  and  (  3. 10)  yields; 


c3    =   -c2X 
cl    =   %/  k 


(3.11) 


Also,  from  Equations  (  3. 6.  a)  and  (  3.10),  we  get: 


v 

(qQ/X  )sinhXdf   -    c2  X  df  sinhXdf 
sinhXdf     +     XdfcoshXdf 


(3.12) 


Now,  for  the  plane  stress  problem,  the  vertical  strain,  €v,  is  related  to 
the  stress  components  by  the  relation: 


dv 


E, 


(  aY    -    vf   oX)     , 


(3.13) 
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where 


Ef      =     modulus  of  elasticity  of  the  foundation 
Ur      =     Poisson's  ratio  of  the  foundation. 


Upon  substitution  of  Equations  (  3.8)  and  (  3.9)  in  (  3.13),  and  integrating 
with  respect  to   y  ,    one  gets: 

cos  Ax      r 

V  (x,y)     -   " —       cl(!  +  vi)  Xsinh  Ay    +   c2(l  +  i/f)A  coshAy 

Ef  L 

+    c3(l  +  ^f)Ay  sinhAy    -    c3(l  -  i/f)cosh  Ay 

+   04(1  +  Vf)AycoshAy   -    c4(l  -  i/f)sinhAy      +  h(x), 

(3.14) 


where 


h(x)     is  the  function  of   x   alone. 


Similarly, 

*x     = 

du 

1 
Ef 

(  ° 
\      X 

■  vi  V 

(3.15) 


where 

u     =     the  deformation  in  the  x-direction 


€x     =     axial  strain  in  the  x-direction. 


Substituting  Equations  (  3.8)  and  (  3.9)  in  (  3.15)  and  integrating  with  respect 
to    x    yields: 

sinAx       r 
u(x,y)    = — ~ c1AcoshAy    +   c2AsinhAy 

+   c3(2  sinh  Ay  +  Ay  cosh  Ay)    +C4(2coshAy 

+   Ay  sinhAy)]    +     g(y)     , 

J  (3.16) 
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where  g(y)  is  the  function  of   y   alone. 

From  Equations  (  3.14)  and  (  3.16),  together  with  the  stress-strain 
relation: 


Qu      +       av     =  TxY       _   2(1  +Vf) 

dy  dx  G  Ef 


yxy     =     "f^    +  ■—  ^^-Txy   ,         (3.17) 


we  get: 

dh(x)     _     _    dg(y) 


dx  dy  C5  <3'18> 


Equation  (  3. 18)  leads  to: 

h(x)     =     C5x    +    C6    ,  (3.19) 

g(y)     =     C5y    +    C?      .  (3.20) 

In  this  particular  problem,  there  is  no  rigid  body  displacement.     This 
implies  that  one  must  have: 

C5     =     C6     =     C?     -     0.  (3.21) 

Substituting  Equations  (  3.11),    (  3.12),  and  (  3.  21)  into  (  3.  14),  and  using 
the  boundary  condition  (  3.6.b)  yields: 

qn  sinh2  A  df 
c2     -   -_ u- 1 .  (3.22) 

X  (Xdf  +  sinh Adf cosh Xdf) 

With  Equation  (  3.22),  one  obtains  from  Equations  (  3.  11)  and  (  3.  12): 

o 
qn  sinh  A  df 

c„     =         — y 1 (3.23) 

**  A(Adf    +   sinhAdj  coshAdf) 
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and 

qQ  sinh  Adf 
X(sinhAdf    +    Adf  cosh  Adf) 


c 


o  (3.24) 

qQdf  sin  h°  Adf 


(sinh  Adf   +    Adf  coshAdf)(Adf   +   sinhAdf  coshAdf) 

Substituting  Equations  (  3.11.b),    (  3.22),    (  3.23)  and  (  3.24)  into  Equation 
(  3.  14)  leads  to  the  following  expression  for  the  vertical  deflection: 

w(x)     =     v(x,  0) 

or, 


w(x)     =     2q   cos  Ax 


or,  in  view  of  Equation  (  3.1): 


o 

sinh    Adf 


Ef  A  (Adf  +  sinhAdf  coshAdf) 


(3.25) 


q(x)     =     %A  (Adf  +  sinhAdf  coshAdf)      .    w(x)  (g  2g) 

2  sinh2  Adf 


If  one  denotes, 

Q(x)     =     2bq(x)  (3.27) 

to  be  the  load  acting  across  the  width    2b    of  the  foundation,  then  one  has 
from  Equation  (  3.26): 

Q(x)     =     k'w(x),  (3.28) 
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where 


b  A  Ef  ( A  df  +  sin  h  Adf  cos  h  Adf) 


sinh2  A  dj 


(3.29) 


The  load  Q,  acting  on  top  of  the  foundation,  can  then  be  expressed  in  terms 
of  the  surface  deflection. 

Consider  next  an  infinitely  long  beam  of  width  2b,  resting  on  top  of  the 
elastic  foundation  subjected  to  the  sinusoidal  load,  p(x),  per  unit  length, 
acting  on  top  of  the  beam  over  the  width  2b,  with  the  foundation  reaction  Q(x) 
defined  in  Equation  (  3.  28)  as  shown  in  Figure  2. 


Elastic  beam 


1 
d 

i 

df 

1 

f^-w^i 


Rigid  Base 


Figure   3#2 
The  load  p(x)  is  expressed  as: 

p(x)     =     PqCosAx 

Also  note  that: 


Elastic 

Foundation 


(3.30) 


Q(x)     =     QqCosAx     =     2bqQcosAx   .  (3.31) 

According  to  the  beam  theory,  the  vertical  deflection  of  the  beam,    w  ,  is 
given  by  the  differential  equation: 
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¥b#     =      P   -   Q-   ,  (3.32) 


where 

E^     =     modulus  of  elasticity  of  the  beam 

Ib      =     moment  of  inertia  of  the  beam  section. 

Substituting  Equations  (  3.28)  and  (  3.30)  in  (  3.32)  yields: 

d4w 
Eb!b   — T    +    k'w     =     Pncos^x     .  (3.33) 

dx^  v 

To  find  the  particular  solution  of  (  3.33),  let: 

w(x)     =     WqCosAx    .  (3.34) 

Substituting  (  3.34)  into  (  3.33),  one  gets: 


w        =     -5 Q .  (3.35) 

A  EbIb    +    k< 


Hence,  upon  substitution  of  Equations  (  3.35)  and  (  3.29)  into  (  3.34), 
we  get: 

2 
p(x)  sinh    A.df 

w(x)    =  — ■ , 

A  EbIb  sinh    Adf   +   bAEr(Adf   +    sinh  Adr  coshAdr) 

(3.36) 

where  p(x)  is  given  in  Equation  (  3.30). 

The  bending  moment  due  to  the  load    p    can  be  evaluated  from  the 

relation: 

d2 
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With  Equations  (  3.36)  and  (  3.37),  it  follows  that: 


M(x)     = ^^ ,  (3.38) 

3  E.pb 


A      +    £r  ^(Adf) 


Eul 


where 


b% 


\df    +    sin  h  Adf  cosh Adf 

<A(Ad  )    =       g •  (3-39) 

sin  h    A  df 


By  means  of  Equation  (  3.38),  it  is  possible  to  calculate  the  bending 
moment  due  to  any  loading  condition  using  the  superposition  principle    and 
the  Fourier  integral.     Recall  that  any  arbitrary  load  p(x)  symmetrical  about 
the  y-axis  can  be  represented  in  the  form  of  a  Fourier  integral: 


/ 


00 


1 

p(x)     =      „•  J  A(A)cosAxdA       ,  (3.40) 

0 


where 

r  oo 


A(A)     =  J  p(v)cos\v  dv     .  (3.41) 

-    00 


Consequently,  the  bending  moment  due  to  any  arbitrary  load  p(x),  expressed 
in  terms  of  the  Fourier  integral,  is  given  by: 

/•   oo 


1 

1 

0 

[A(  A)cos  Ax  ]    AdA 

TT 

v3        Efb     , 
EbTb 
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B.      An  Infinite  Beam  Subjected  to  a  Uniform  Strip  Load 

Suppose  that  the  beam  is  subjected  to  a  uniform  strip  load  of  intensity 
p    and  width  2  €  as  shown  in  Figure  3. 


Elastic  Beam 
x 


Elastic  Foundation 


Rigid  Base 


In  this  case,  one  has: 


Figure     3.3 


p(x)     = 


p  , - €<X<€ 


0       elsewhere 


(3.43) 


Hence: 


A(A) 


/ 

-  e 


p  cos  Av  dv 


or, 


A(A)     =     2p€ 


sin  A€ 
A€ 


(3.44) 


Substituting  (  3.  44)  in  (  3.  42),  and  letting 

1/3 

r  Eb!bi 

a 


[Ffb     J 
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(3.45) 


and 

a  =     aX  ,  (3.46) 

where      a     =     a  fundamental  length 

a     =     a  dimensionless  parameter, 

leads  to: 

■00 


/oo 
sin(a< 
a3  + 


2pa  /       sin(a£/a)cos(ax/a) 

M(x)  — /      $ ;    ;       ; — —L      da      .  (3.47) 

QJ         ad  +  (/r(adf/a) 


Similarly,  it  can  be  shown  that  the  corresponding  expressions  for  the  deflec- 
tion w(x),  the  shearing  force  V(x)  in  the  beam,  and  the  contact  pressure  Q(x) 
(load  per  unit  length)  at  the  interface  of  the  beam  and  the  foundation  are: 


r 

/        sin  (a  €/, 


w(x)     =    f^  J        sin(a€/a)cos(ax/a) 


77  J 


00 


V(x)  -^       J     Bin(a€/a)sln(ax/a)       .    dfl  (3>49) 

0J  a  6  +  ijj  (odf/a) 


/oo 
^(adf/a)sin(a  t/a)cos(qx/a)      ,  dfl 
a  fa3    +  «£  (adf/a)l 


00 

Q(x) 


C.      An  Infinite  Beam  Subjected  to  a  Concentrated  Load 

The  elastic  solution  to  the  problem  of  an  infinite  beam  subjected 
to  a  concentrated  load,  P,  acting  along  the  y-axis  can  be  obtained  as  a 
limiting  case  of  the  uniform  strip  load  problem.      As  the  loaded  region 
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-0,  one  has: 


2p£       , 


lim  sin(a€/a) 

€ ►  0        a  €  /a 


Consequently,  the  expressions  for  w(x),  M(x),   V(x),  and  Q(x)  due  to  the 
load  P  take  the  forms: 

3  "°° 


cos  (ax/a) 

— y—T —J — •    da  (3.51) 

a  fa    +  <A  (adf/a)  | 


/oo 
a  cos  (ax/a) 
a3  +  ^(adf/a)  da  <3-52> 

/OO            g 
a     sin  (ax/a) .    da 
„3      +i//    f«rWfl\ 


0 


a°    +  <£  (adf/a) 


(3.53) 


Q,x')     =    -£-       /"»(«0/a)°0B<ax/.)       .    ,„ 

0*7  a      +i/a  (adf/a) 


D.      An  Infinite  Beam  On  A  Semi-Infinite  Foundation 

This  is  a  limiting  case  of  an  infinite  beam  resting  on  a  finite  depth 
foundation  solutions.     Rewrite  Equation  (  3.39)  in  terms  of  the  new  para- 
meter a  : 

<Madf/a)     =    ~ g  df/L+  1^  (°df/a)  C°Sh  (gdf/a)        •       (3-39a) 


sin  h     (adf/a) 
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It  is  seen  that  as  cL*  »    oo  ,  \\j  (adf/a)  >   1  as  a  limit.      Hence, 

for  the  case  of  a  semi -infinite  foundation,  Equations  (  3.47)  through  (  3.54) 
are  still  valid  upon  substitution  of  ip (adf/a)     =     1.0. 

It  is  noted  that  the  expressions  thus  obtained  for  the  deflections, 
bending  moments,  shearing  forces,  and  the  contact  pressures  are  very 
similar  in  character  to  those  obtained  by  Biot(l),  Vesic  (8),  Drapkin  (2), 
and  Selvadurai  (5).      The  solutions  for  the  case  of  plane  stress  problem 
of  an  infinite  beam  on  a  semi-infinite  foundation  are  exactly  the  same  as 
those  obtained  by  Biot. 

For  a  plane  strain  problem,  one  simply  replaces    Ef  with  Ef/(1  -  ut  ) 
in  the  expression  for  the  fundamental  length  'a'  as  defined  in  Equation  (  3.  45). 
That  is,  for  a  plane  strain  problem,  one  has: 

2,  „  T  -|l/3 


a     = 


(1  -  vf")  EbIb 


(3.55) 


Efb 

With  'a'  defined  in  Equation  (  3.55),  the  expressions  for  the  deflections, 
bending  moments,  etc.  for  the  plane  stress  problem  can  then  be  used  for 
the  plane  strain  problem. 

E.      Numerical   Evaluation  of  Related  Integrals 

A  numerical  technique  that  can  be  used  to  evaluate  the  finite  integrals 
appearing  in  the  preceding  section  is  shown  here  in  detail  for  the  case  of  the 
beam  subjected  to  a  concentrated  line  load.     The  same  technique  also  applies 
to  the  other  cases. 

Consider,  for  example,  the  expression  for  the  bending  moment  of  the 
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beam  due  to  a  concentrated  line  load  given  by  Equation  (  3.52).     Transform 
the  equation  to  the  form: 


M(x)     =     _1_ 
Pa  it 

0 


/oo 


ea  a  cos  (ax/a) 


a3    +  if/  (adf/a) 


da  (3. 56. a) 


or 


00 


T*a       =     -jj-        }        e   a  f  (a,  x/a,  df/a)  da  (3.56.b) 

0 

In  this  form,  the  integral  can  be  evaluated  by  a  Gauss-Laguerre  quadrature 
formula . 

A  computer  program  was  written  to  evaluate  Equation  (  3.  56)  for 
various  values  of   x/a    and    df/a,  using  the  32-point  Gauss-Laguerre  quadra- 
ture formula.    As  a  check  for  accuracy,  one  might  use  Equation  (  3.  52)  to 
evaluate  the  maximum  moment  (at  x  =  0)  when   df  »    oo  ,  which  results 

in: 

00 


Pa  / 

IT  J 


^r  ;  -sir  da  <3-57) 

0  a  L 


This  integral  can  be  integrated  exactly  to  get: 

Mmax     =     °-385  Pa  (3.58) 

as  compared  to  the  computer  output: 

Mmax     =     0-382  Pa     .  (3.59) 
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Hence,  the  error  is  in  the  order  of  1%  or  less. 

Figure  4  shows  a  series  of  normalized  bending  moment  curves  ob- 
tained from  numerical  evaluations  of  Equation  (  3.  56). 

From  Equation  (  3.  56),  the  maximum  moments  for  different  values 
of  df/a  are  given  by: 


W 


00 


Mmax  1  /  a    dct  (3.60) 

0  a      +    ip  (adf/a) 


If  the  numerical  results  of  Equation  (  3.  60)  are  plotted  on  logarith- 
mic paper,  it  is  found  that  for  0.1  l_  df/a  l  1.5,  Equation  (  3.  60)  can  be 
represented,  with  an  approximation  on  the  order  of  0.  5%,  as: 

0.2438 


M  /df 

max I 


=     0.2904  1 1  ,     0.1    4   df/a   <    1.5 


Pa  \  a 


or 


df 


0.2438 


Mmax     =     °-2904Ui  Pa  (3'61^ 


Now,  from  the  fundamental  theory  based  on  Winkler's  assumption, 

it  is  given  that: 

/      i  \l/4 
Mmax     =     0.353554  P  (-^)  (3. 62.  a) 


k     / 


where  k     =     modulus  of  subgrade  reaction    . 

Hence,  in  order  to  obtain  the  same  maximum  bending  moment  by  the 
approximate  theory  as  by  the  exact  one,  one  might  have  to  use  a  modified 
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value  of   k   found  by  equating  Equations  (  3.  61)  and  (  3.  62.  a),  i.e 

1/4 


0.353554  P 


^) 


v  0.2438  .  .   1/3 

df  \  /e*~ 

[1T  P\E, 


Vb\ 
Sfb    y 


or 


Ef 


=     2.197 


df 

a 


-0.9752 


0.1   4   df/a  4   1.5 


(3.63) 


where    a    is  given  by  Equation  (  3.45). 


It  should  be  noted  that  Equation  (  3.  62.  a)  is  derived  from  the  general  expres- 
sion for  the  bending  moment: 


M(x)     = 


-qrg-    e  (cos    fix    -    sin  fix) 


(3.64) 


from  which 


M 


max 


40 


(3.62b) 


where 


0 


4  EvJ 


bxb 


1/4 


So,  in  view  of  Equations  (  3.61)  and  (  3.62.b),  one  gets: 


fia     =     0.8609 


df. 

a 


-0.2438 


where 


0.1   4   df/a    <    1.5 


(3.65) 


(3.66) 


Hence,  knowing  fi,  one  can  also  compute    k   from  Equation  (  3.65),  which 


gives : 
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k     =     4  04EbXb  (3.67) 


Equations  (  3.  66)  and  {  3.  67)  are  suitable  forms  for  evaluating  fi  and    k 
respectively.      Equation  (  3.  66)  has  been  evaluated  numerically  and  is 
plotted  in  Figure  5. 

Knowing  the  stresses  and  the  contact  pressure  in  the  uncracked 
beam,  it  is  then  possible  to  use  these  stresses  as  known  boundary  condi- 
tions to  the  approximation  of  the  stresses  in  the  beam  having  a  centrally 
located  vertical  crack . 
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IE.      AN  ANALYSIS   OF  STRESSES  AND  THE  OPENING -MODE 
STRESS  INTENSITY      FACTOR    FOR  A    CRACKED  BEAM 

A  method  of  analysis  consists  of  finding  a  stress  function,  X  ,  that 

4 
satisfies  the  biharmonic  equation    V     X    =     0,  and  also  the  boundary  con- 
ditions along  the  boundary  of  the  cracked  beam.     The  biharmonic  equation 
and  the  boundary  condition  along  the  crack  are  satisfied  by  the  Williams 
stress  function  (9),  which  in  the  case  of  symmetry  of  loading,  is  given  by: 

00  / 

X(r,6)     =      X  \ M)11"1  d2n-l    r(n+l/2)    [  -cos(n-3/2)  0 

n=l,2,...     ( 


2n-3  ,„      t 

+    "2nTT     C0S  (n  +  l/2)(9] 


.(n+1) 


+    (-l)n   d2    r[       }   [-cos(n-l)^  +   cos(n+l)(9] 


(3.68) 


P/2  6D 


Crack 


T77TTT7TTTTT7^JT7T7T7T71TT7T7 


Rigid  Base 


Elastic  Beam 


Elastic  Solid 
Foundation 


Figure  3.  6      Cracked  Beam  on  Elastic  Foundation 
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According  to  the  coordinate  system  shown,  the  stresses  in  terms  of  X  are: 

d2X  d2X  2.  d2X         sin0  cos  $ 

°v     "        ^    2     _        „    2      cos  ^  "    2 


y  9x"  ar  30dr  r 


gy  sin2 e  ,      o      ax  sine  cos  g 

+      — +       Z     — — 9 

ar  r                          dd                  * 

+     g  y  sin2^ 

a«2  r2 


o-       =     Jg^.a    -ggx_    sin2g+2_a^     sine  cos ^AKcos2^ 


av*        ar*  aear 


^2  2„ 

9    a  X         sin  e  cos  e     .       o   X       cos  0 

-    &     r  - —       T ■+-    — r — tj-     rr 

r 


do         r^"  ae2 


a2X  .    n        a         d2X  cos  20  d2X 

xy  ax  ay  5rz  r  O0dr 


+       g  X1       sing  cose         +     d  X  sin  flcose 


ae2        r2  a 


r  r 


+ 

dx 

ae 

cos  2e 

r2 

where 

ex 

00 

V 

J,.„ 

(3.69) 


ar 


X  ]  (-l)n_1  (n+1/2)   d2n_i    r(n-l/2)     [-cos(n-3/2)e 

n=l,2,  .  .  .  ! 

2n-3 
+  ^nTf   C0S(n+V2)e]   +    (-l)n  (n+l)d2n    rn 

[-    cos(n-l)e     +   cos(n+l)e]  \  (3.70. a) 
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a2v  °°  I 

~d~^    =       £  |(-l)n-1  (n2-l/4)d2n_1r(n-3/2)     [-cos(n-3/2)0 

n=l,2,  ..  .1 

+   lv+f     cos(n+l/2)(9]  +    (-l)n  (n+l)d2nr11-1  (3.70.b) 

[-  cos(n-l)0    +    cos(n+l)#] 

dX 


e$ 


^2  I  ("if"1  d2n-l  r(n+1^2)    [(n-3/2)  sin(n-3/2)0 

n=l,2}...    ( 


-    (n-3/2)  sin  (n+l/2)  fl]  +    (-l)n  d2n  rn+1  (3.70.C) 


[(n-1)  sin  (n-l)0     -    (n+1)  sin  (n+l)fl]  > 


,2. 


oo 


-^4"    =    E  {("I)11"1  d2n-l  r(n+l/2)    [n-3/2)2Cos(n-3/2)0 

£0  i£f,2,...    ( 


-    (n-3/2)  (n+l/2)  cos  (n+l/2)0]  (3.70.d) 


+  (-l)n  d2n  rn+1    [(n-l)2cos(n-l)0     -    (n+l)2cos(n+l)(9] 


00 

d2x 


eee 


=      J2  \  (-1)11"1  d2n-l  <n+1>  r^"1/2)  [(n-3/2)  sin  (n-3/2  )0 

n=l,2, ...    ( 

-    (n-3/2)  sin  (n+l/2)fl]      +     (-l)n  (n+1)  d2n  r11  (3.70.e) 

[(n-1)  sin(n-l)0   -    (n+1)  sin  (n+l)0]> 


If  it  is  assumed  that  the  presence  of  the  crack  does  not  have  any  signi- 
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ficant  effect  on  the  stresses  beyond  the  boundary  BC,  as  well  as  along  the 
interface  of  the  beam  and  the  foundation,  the  stresses  along  those  boun- 
daries can  then  be  closely  approximated  from  the  conditions  of  the  uncracked 
beam.     Taking  symmetry  into  consideration  and  considering  the  portion  of 
the  beam  cut  along  BC  as  a  free  body,  one  has  the  following  additional 
boundary  conditions: 


Figure  3.7 


P/2eB 


Crack 


Along  AB: 


T         -     0 

xy 


X 


7T6B 


00 

/- 


i/r(adf/a)  sin(flt£/a)  cos(ay/a) 
0  a[a3    +  ,/,(adf/a)] 


da 


Along  BC: 


xy 


°y 


=    r 


xy 


/oo 
a  sin(q  e/a)  sin(a  //a)     $ 
a'6    +   il/(adf/a) 


(//■(adf/a) 


=     <7' 


-/x\     (pa  )       /     sin(q  e/a)  cos(q  //a)     da 


0 
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(3.71) 


(3.72) 


where 


f  (x)     =     ~ T3-     cd  +  xd  -  2xc  -  c2  -  x2 


i    r        d 
g(x)  =   "lb~LX  "  2  +  c 


]■ 


Along  CD: 

T 

xy 

=     0 

a 

X 

=     0 

Along  DE: 

T 
xy 

=     0 

°x 

=     -  P/2 €  B 

(3.73) 


(3.74) 


The  problem  will  be  explicitly  solved  if  we  can  find  the  coefficients  of  the 
Williams  stress  function,  dj  ,  i  =  0,1,2, ...  ,  such  that  Equation  (  3.69) 
satisfies  the  boundary  conditions  given  by  Equations  (  3.  71)  through  (  3.  74). 

As  for  the  opening-mode  stress  intensity  factor,  K^,  we  have,  accor- 
ding to  Irwin  (3),  in  the  immediate  vicinity  of  the  crack  tip  (as  r  approaches 
zero): 


K 


V27T: 


cos     -|-    (1  +  sin  -@-   sin  -g-  ) 


(3.75) 


while  in  terms  of  the  Williams  stress  function, 

cos   JL  (1  +  sin  -f-   sin   ^~ 


V^~ 


(3.76) 


hence, 


Kt      =      -     V2~7T 


1 
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(3.77) 


Therefore,  only  the  value  of  the  first  coefficient,  d^  f  is  needed  to  determine 
the  value  of  the  stress-intensity  factor,  Kj  . 

A.      An  Approximation  of  the  Stress-Intensity  Factor  by  the 
Boundary  Collocation  Procedure 

Gross  and  Srawley  (6,  7)   used  a  boundary  collocation  procedure  in  con- 
junction with  the  Williams  stress  function  in  the  approximation  of  stress-intensity 
factors  for  single-edge-notch  beams  in  bending  and  combined  bending  and  ten- 
sion.    Later,  Majidzadeh  and  Ramsamooj  (4)   adopted  the  same  method  to 
approximate  the  stress-intensity  factor,  Ki ,  for  beams  having  a  centrally 
located  crack,  supported  on  a  three-dimensional  semi-infinite  elastic  solid, 
and  subjected  to  a  centrally  located,  vertically  concentrated  load. 

The  same  technique  is  employed  here  for  solving  the  two-dimensional 
problem  of  an  infinitely  long  beam  having  a  centrally-located  vertical  crack 
at  the  bottom,  supported  on  an  elastic  solid  of  any  thickness,  and  subjected  to 
a  uniform  strip  load  (a  concentrated  line  load  is  merely  a  limiting  case). 

The  boundary -collocation  procedure  consists  of  solving   2m    simul- 
taneous algebraic  equations  corresponding  to  the  known  values  of  stresses 
and  displacements  at  a  finite  number  of  stations,    m  ,  along  the  boundary 
A  B  C  D  E  (Figure  7)  of  the  free  body  of  a  cracked  beam.     The  solution  yields 
the  values  of  the  first    2m    coefficients  of  the  Williams  stress  function.     For 
the  present  purpose,  only  the  first  coefficient,  d-, ,  is  needed  in  the  determina- 
tion of  the  stress-intensity  factor,  Kj  . 
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For  given  values  of  geometrical  dimensions  and  loading  conditions, 
the  value  of  d-,  varies  somewhat  with  the  number  of  terms,  2n  ,  of  the  stress 
function,  the  number  of  boundary  stations,  m,  and  the  length    /  .      It  has  been 
found  that  for  a  given  2n,  the  solution  for  d-^  converged  rapidly  with   l/d  =  1.0 
and  with  the  stations  spacing  within  the  range  d/8  to  d/16.     To  further  opti- 
mize the  solution,  it  has  been  recommended  (4)  to  use  the  linear  least-squares 
optimization  scheme  in  solving  the  problem;  using  more  equations  than  un- 
knowns.   Using  this  particular  technique,  it  was  noted  that  the  optimal  num- 
ber of  dj  terms,  2n,  varies  somewhat  with  each  particular  beam's  dimensions 
as  well  as  with  the  crack  length  over  beam  depth  ratio;    however,  their  effects 
on  the  values  of  K  j  is  practically  small. 

A  value  of  Kj  obtained  by  this  method  corresponds  to  a  particular  set 
of  values  of  P,  d,  / ,  c,  m,  and  n.     For  application,  the  results  are  better 
expressed  more  generally  in  terms  of  the  dimensionless  quantity, 

Kl'    =     V^maxVd")  (3.78) 

where  0"max    is  the  maximum  tensile  stress  at  the  bottom  of  the  uncracked 
beam  on  elastic  solid.    With    l/d    >  1.0,  it  was  found  that  K^'  is  practically 
dependent  on  the  c/d  ratio  for  a  set  of  m  and  n. 

For  computation  convenience,  a  beam  section  and  loading  conditions 
shown  below  was  selected  for  the  approximation    of  Kj '  for  various  c/d 
ratios: 


158 


of 


p  =10  psi 


r 


€  =  0.5 


"symmetry 


/  =  2< 


EE 


Spacing  between  each  station  = 

1/8" 


d  =  2' 


Figure  3.  8 

It  was  found  that  the  optimum  number  of  d^  terms,  2n  ,  for  various  c/d 
ratios,  lies  between  60  and  80.    Hence,  for  each  particular  c/d  ratio,  the 
value  of  Kt  '  was  obtained  by  taking  the  average  of  three  values  computed 
for  2n  equal  to  60,  70,  and  80. 

The  results  of  the  computations  are  plotted  in  Figure  9.     To  verify 
the  validity  of  this  normalized  relation,  different  beam  sizes  were  selected 
to  compute  Kj '  for  a  particular  c/d  ratio.     The  results  agreed  very  well. 
The  computer  program  used  in  the  boundary  collocation  method  is  presented 
in  Appendix  A . 

B.      Applications  to  Analysis  of  Fracture  Propagation 

For  a  particular  c/d  ratio,  the  value  of  crmax  under  a  set  of  geometri- 
cal and  loading  conditions  must  be  known  in  order  to  compute  the  absolute 
value  of  Kj  from  the  normalized  Kj'.     For  a  rectangular  cross  section,  one 
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Figure  3.  9      Relation  Between  the  Normalized  Stress-Intensity  Factor 
and  the  Relative  Crack  Length. 
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has: 


6Mmax 
'max     =     "US-  <3.79) 


where  M  is  the  maximum  bending  moment  in  the  beam.    M^,ov  can  either 

max  &  max 

be  computed  from  Equation  (  3.47)  for  a  uniform  strip  load,  or  from  Equation 
(  3.52)  for  a  concentrated  load.    However,  this  direct  method  for  evaluating 
M  requires  a  numerical  integration  which  is  only  feasible  with  the  aid  of 

a  computer.    A  practical  approach  in  the  evaluation  of  arnav  is  proposed  in 

XXX  dL  a 

the  following. 

It  has  been  shown  that  the  bending  moment  in  a  beam  resting  on  an 
elastic  solid  of  a  finite  depth  and  subjected  to  a  concentrated  load  may  be  com- 
puted from  the  conventional  beam-on-elastic-foundation  formula;  i.e.  : 


M(y)    -    — -     e   0y  (cos  0y    -    sin  0y)  (3.64) 


if  fi  is  taken  as: 

-0.2438 


a    =    0.8609 


df 


,   0.1<  df/a  4  1.5    .  (3.66) 


Consequently,  one  gets: 


Mmax     "    T7  <3-62'b> 

Substituting  Equations  (  3.  66)  and  (  3.62.b)  into  Equation  (  3.79)  yields: 

rdfi( 

1.7424  1—  ,   0.1  4  df/a  <  1.5    .        (3.80) 


a  Bd2  rd,1°-2438 

p  max 


Pa 
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It  is  noted  that  Equation  (  3.80)  is  applicable  for  0. 1  <  df/a  <  1.  5,  which  is 
considered  to  be  within  the  practical  range  of  df/a  in  an  experimental  set-up. 
For  a  special  case  of  semi-infinite  foundation,  it  can  be  shown  that: 


a  Bd2 

p  max 

=     2.3094  (3.81) 


Pa 

The  plot  of  Equation  (  3.  80)  is  presented  in  Figure  10. 

Thus  far,  one  has  been  able  to  approximate  the  maximum  tensile  stress 
0"p  max  m  the  beam  due  to  a  concentrated  load,  P.     To  extend  the  application 

to  a  uniform  loading  conditions,  a  set  of  numerical  integrations  of  M    „     from 

°  °  max 

Equation  (  3.  37)  were  carried  out  for  0    4  e/a  4  2.0.     The  results  were  ex- 
pressed in  terms  of    (7S  max    /  0"p  max  »  wnere   O 's  max  is  the  maximum  ten- 
sile stress  due  to  a  uniform  load  and  O     max  is  the  maximum  tensile  stress 
due  to  a  concentrated  load.      Figure  11  shows  the  plot  of  as  max/cr     max 
against  e/a  ratios."  For  a  particular  e/a  ratio,  the  corresponding  ratio  of 
<7G  v>tqV/  (7„  ™o^  can  then  be  read  from  Figure  9.    Knowing  this  ratio  and 

o    IIlcLA  p   IllctX 

an  „„._  from  Equation  (  3.  80)  or  Equation  (  3.  81),  or  from  Figure  10,  one 
can  then  compute    as  max. 

Having  obtained  the  maximum  tensile  stress  for  the  unc racked  beam, 
the  absolute  value  of  the  stress-intensity  factor,  Kj,  may  then  be  computed 
from  Figure  9.      However,  instead  of  directly  reading  the  value  of  Kj  '  for  a 
c/d  ratio  off  of  Figure  9,  the  following  empirical  equation  was  found  to  be  a 
compact  expression  of  the  present  results  for  0.14  c/d  <  0.  6, 

Kj'      =     6.898(c/d)    -    17.425(c/d)2    +    22.438(c/d)3  (3.82) 
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Equation  (  3.  82)    is  plotted  in  Figure  3.12.    K  may  then  be  approximated 
from  Equation  (  3.82). 
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KT    2 


^^■^^^^ 


X   Collocation  results 


KT= 


*L 


1      o     d 

"max  V* 


1/2 


=  6.898  c/d-17.425[c/d]  +22.438[c/dJ 


.2  .3  4 

Relative  crack  length,  c/d 


•5 


Figure  3. 12        Approximation  of  the  Normalized  Stress-Intensity  Factor 
as  a  Function  of  the  Relative  Crack  Length 
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IV.        GEOMETRICAL  DESIGN  OF  BEAM  SAMPLES 
AND  EVALUATION  OF  Kj 


A.  A  Criterion  For  the  Length  of  a  Beam  Sample 

It  is  known  from  the  approximate  theory  based  on  Winkler's 
assumption  that  if 

0  I    >     5      ,  (3.83. a) 

where  /     =     the  length  of  the  beam, 

then  the  end  conditioning  forces  at  either  end  of  the  beam  of  finite  length  / 
will  have  an  insignificant  influence  at  the  other  end  and,  therefore,  all  the 
formulas  derived  from  the  infinitely  long  beam  can  be  applied  to  this  case. 
In  this  case,  knowing  the  value  of  fi  from  Equation  (  3.66),  the  minimum 
length  of  the  beam  sample  then  should  not  be  less  than: 

/     =     5/fi  (3.83.b) 

B.  Criteria  For  the  Cross  Section  of  the  Beam 


According  to  the  ASTM  Specification  C31-69,  it  is  recommended 
that  for  flexure  tests,  a  beam  should  have  the  minimum  dimension  of  three 
times  the  maximum  nominal  size  of  the  coarse  aggregate  in  the  mix.     Further- 
more, the  depth  of  the  beam  should  not  exceed  four  times  its  width.    Hence, 
the  beam  cross  section  should  be  such  that: 

B     =     2b   >  3  (max.  size  of  aggregate)  (3.84) 
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and 

B  <    d    <  4B  (3.85) 

where  B     =     the  width  of  the  beam 

d      =     the  depth  of  the  beam.  . 

C.      A     Criterion.  For  the  Determination  of  the  Load  Magnitude 

In  order  to  prevent  excessive  rutting  of  the  beam  top  surface  under 
the  load,  as  well  as  to  assure  that  the  maximum  tensile  stress  at  the  bottom 
of  the  beam  surface  is  below  the  yield  stress,  it  is  recommended  that  the 
applied  load  should  be  of  such  a  magnitude  as  to  cause: 


Vx    £     °'5     °y  <3"86> 


where  cr         -     yield  stress  in  tension. 


Knowing  the  beam  and  the  foundation  dimensions,  together  with  the 
related  material  properties,  one  can  then  select  the  required  magnitude  of 
the  applied  load  from  Figures  11  and  10  respectively. 

D.      Procedure  For  the  Design  of  a  Beam  Sample 
a)    Selection  of  the  beam  dimension: 
1.       Evaluate  the  values  of 

E^     =     modulus  of  elasticity  of  the  beam 

Ef      =     modulus  of  elasticity  of  the  foundation  from  the  experiments 
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2.  Select  a  minimum  width  of  the  beam  from 

B     =     3  (maximum  nominal  size  of  aggregate) 

3.  Select  a  series  of  beam  depths,  d,  such  that 
B  4  d  4  4B 

1  O 

4.  Find  the  section  property   Ib     =    YI  Bd°  and  from  which 
calculate: 

a  =  ( -ETb-y      • where  b  =  B/2 

5.  Select  a  foundation  depth  d^  and  then  compute  df/a.     Read 
off  the  curve  in  Figure  5  the  corresponding  value  of  0  a. 
Then  compute  0  . 

6.  The  required  length  of  the  beam  is  computed  from: 

/    =     5/0 

7.  For  optimum  solution  for  each  particular  E^  and  Ep  try 
different  beam  depths  as  well  as  different  foundation  depths, 


Example ;         Given  the  following  properties 


Beam  size:     3"   x    3"  Ib    =  A  (3)4    =    6. 75  in. 4 


Eb     =     65,000  lb. /in.2,      Ef     =     250  lb. /in.2 
/EkIkW3        /  65,  OOOx  6.75  \  !/3 

a  =  Knr      ■-(  250  x  1.5         =  10-54in- 
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For   df/a    =    1.5" 


df/a     =     1.5/10.54     =     0. 


142 


From  curve,  one  gets  with  df/a    =    0. 142, 
0  a     =     1.375 

8  =     1.375/10.54    =    0.1305  in.-1 


/   =     5/0      =     5/0.1305    =    38 


in. 


Similarly,  for  df  =  2",  one  gets  /  =  41" 
df  =  3",  one  gets  /  =  45" 
d      =    4",   one  gets     /  =    48.5" 


8  5CH 

o 

'"  45  - 


f  40 

CD 


a   35  - 


30 


i 

/ 


Eb=  65,000   lb/  in: 


Ef  = 


/ 


250   lb./ in 

Beam    size :     3"x  3" 


2 


5 


~r 
6 


-»^  df  in. 


12         3        4 

Figure  13.    Relation  between  the  required  beam  length  and  foundation  depth 
b)     Selection  of  the  magnitude  of  the  applied  load 

1.  Select  the  length,  2  e  ,  of  the  cushion  used  to  transfer 
the  load  from  the  loading  piston  onto  the  beam  surface 
(the  width  of  which  should  equal  the  width  of  the  beam) 
as  shown  in  Figure  14. 
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Cushion 


b,b 


^ 


Section  A-A 


Cushion 


Beam 
Sample 


T7C 


Foundation 


Rigid    Base 


FIGURE    14 


Compute   e/a  and  use  this  value  to  read  off  the  curve  in  Figure 
11  the  magnitude  of  ag  max  in  terms  of  ap  max,  denote  this  by: 

ap  max  as  max  al  ap  max' 

when  a^is  the  value  read  off  from  the  curve. 

3.        Select  a  desired  value  of  ap  max  such  that, 


ap  max  a2  ax 


a2<  1 


where         <7y       =     yield  stress  in  tension 


a  desired  load  factor 


Compute  the  value  of  a'    „,„,,»  i-e- 

L)    III  d-A. 


p  max 


al°2    °\ 


5.        Compute  the  value  of  df/a  and  read  off  from  Figure  10 
the  corresponding  value  of  a    tviqv  Bd  /Pa,  say 

ap  max  Bd 

=       a 

Pa  3 
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6.        Substitute  the  value  of  a'p  max  f°r  ap  max  in  tne  above  expression. 
The  required  magnitude  of  the  applied  load,  P,  can  then  be  com- 
puted from: 


Example ; 


i^'  ";\     """"" 


a^ao  a      Bd 


a3    a 

7.  Check  for  the  contact  pressure  between  the  cushion  and  the 
beam  surface  from: 

p     =     P/2CB, 
where  B  is  the  width  of  the  cushion.     The  pressure  p  should  be 
considerably  less  than    cr      in  order  to  prevent  excessive  rut- 
ting;   it  is  recommended  that 

p    <    0.5  CTy 

8.  If  the  contact  pressure  p  is  considerably  large,  one  can  either 
select  a  new  dimension  of  the  cushion  and  repeat  the  entire  cal- 
culation, or  simply  select  a  new  load  factor,   a2>  and  repeat  the 
computation  from  step  6. 


Beam  size:     B    =    3"  ,    d    =    3" 

Eb     =     65,  000  lb/in?      ,      Ef     =    250    lb/in. 2 
2  €     =     2  in. 


df       =     1. 5  in. 

a       =     150  lb/in.  2 
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One  has  from  the  preceding  example: 
a        =     10.  54  in. 

df/a    =     0. 142 

e/a      =     0.095 
With  e/a    =    0.  095,  we  get  from  Figure  11: 


tf'p  max         as  max         0"  95  ap  max 


Select    a    max     =     0. 5    a       ,  i.e.  a2    =    0.5, 
one  gets: 

a'p  MX     =     0.95    x    0.5    x    150     =     71. 5  lb/in 

With  df/a    =    0. 142  one  obtains  from  Figure  10: 

°p  max    -° 

—     =     a3       =     1.09 


:~  2 


Pa 


Hence,  substitution  of  a'     ^,ov    for    a  in  the  above 

'  "  p  max  p  max 


equation  leads  to: 


71.5    x    33 


P     = 1LJ1 =     168    lb. 

1.09    x    10.54 


Next,  check  the  contact  pressure,  assuming  that  the  cushion 
has  the  same  width,    B     =     3"    , 


26B  2x3 


!§%-     =     28    lb. /in.2 


«  <7y      o.k. 
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c)        Determination  of  Kt  value  for  a  particular  crack  length: 

1.  Compute  for  each  particular  crack  length  c 
the  ratio  c/d. 

2.  For  each  c/d  ratio,  read  off  from  Figure  12  the 
corresponding  Kj'  value. 

3.  Let    (Tni!,Y=    o~'  a-,  a0  0%,  as  obtained 

max       «  p  max  1     2      y 

from  the  preceeding  example. 

4.  Compute  for  each  crack  length  the  corresponding 
value  of  Kt  from  the  following  equation: 

Ki     =        a-max    '  Vd~     K[ 

Example:         Using  the  same  data  as  given  in  the  preceeding  example,  sup- 
posing that  the  value  of  Kj  for  a  crack  length  c  =  1.  05  inches  is  required. 

1.  c/d     =     1.05/3.0    =    0.35. 

2.  With  c/d  =  0.35,  from  Figure  12  the  corresponding 

Kj'  =  1.25 

3.  From  the  preceeding  example, 

°*max     =    «T!  =    71. 5  lb/in2 

max  «  p  max 

4.  With  d  =  3.  0  in. , 

KI     =      %ax     V^~       Kjt 

=     71.5    VY        (1.25) 

-3/2 
=     154.  8  lb-in. 
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APPENDIX    A 


Description  and  Use  of  the  Collocation  Computer 
Program  for  the  Calculation  of  the  Stress-Intensity 
Factor,  K-, ,  for  Two  Dimensional  Problems  of  a 
Beam  With  a  Centrally- Located  Finite  Crack,  Res- 
ting on  an  Elastic  Solid. 
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IA.      GENERAL  DESCRIPTION 

This  program  uses  the  Williams  stress  function,  X ,  in  conjunction 
with  the  boundary  stresses  from  the  elastic  theory  of  beams  on  elastic  founda- 
tion to  compute  an  opening-mode  stress-intensity  factor,  K^,     and  then  nor- 
malizes it  to  K-l'.      The  method  chosen  to  solve  the  problem  is  the  "boundary 
collocation  of  stresses"  ,  which  consists  of  solving  a  set  of  2m  simultaneous 
algebraic  equations  corresponding  to  the  known  values  of  stresses  at  a  finite 
number  of  boundary  stations    m    of  the  cracked  beam. 

Since  it  is  impossible  to  deal  with  the  exact  form  of  the  Williams 
stress  function  given  by  Equation  (       46)  in  the  text  (involving  an  infinite  sum 
of  series),    only  the  first    2n   terms  of  it  were  included  in  the  analysis,  i.e.  , 
the  Williams  stress  function  was  approximated  by. 

NOOFM     (  fn+i/9*      T 

*    =        Z  (-D         d2n-l  r(       7   }       -cos(n-3/2)0 

^     n=l  (  L 

2n-3  ,       1 

+  ~2nTf       cos(n  + 1/2)6)J  (A-l) 

+  (-l)nd2nr(n+1)    |-cos(n-l)0    +   cos(n+l)0|  > 

The  corresponding  stress  components  are  given  by: 

a     =      32*     sin2e    _    2    &2*        sin  e  cos  0        +       d*      cos   0 
x  a  r2  dOdr  r  dr  r 

2dX       sing  cos  0        ,      dV         cos2  d 

a      =  -^   cos2d_2^h_     sing  cos  8      +    6X__    sin2^  (A 

y  dvl  dddr  r  dr        r 

2dX        sine  cosfl  j.       d2X~      sin20 


d$  r2  d&2  r2 


178 


T       =   -  JUL-  sin  e  cose  -  -At    cos  20       +    A_     sing  cos  0 
xy  drl  66dr  r  ae  r 


+      gx      sin  0  cos  0 +      ex       cos  2  0 

dr  r  a«  r2 


where 


-  NOOFM  f 


n-1,2,  ..  . 


2-  NOOFM 

ar2  2-, 

n=l,2... 


-cos(n-  3/2)0       +  2n+1      cos  (n  +  1/2)0 

+  (-l)n   d2n    r11      |- cos (n- 1)*    +  cos(n+l)e~|  |  (A  -3. a) 

J  (-l)^1)  (n2  -  1/4)  d2n_1  r(n-3/2)     [-  cos  (n-3/2)  9 


2n-3 


2n+l 


cos(n+l/2)0|    +    (-l)n    (n+l)nd2nr(n  ^  (A  -  3.b) 


-cos(n-l)0    +   cos(n+l)$ 


_  NOOFM 

dX 


69 

n=l,2, 


£  |    (-l)(n_1)  d2n_1  r(n+1/2)    J  (n-3/2)  sin  (n-3/2  )* 


-  (n  -  3/2)  sin(n+l/2)0  |   +    (-l)n  d2n  r(n+1)  (A  -  3.c) 

(n-l)sin(n-l)0    -    (n+l)sin(n+l)0  |  > 

_  NOOFM     / 

^f     =        E  (-l)(n_1)  d2n_!  r(n+l/2)     r(n-3/2)2cos(n-3/2)0 

n=l,2,...l  L 

-  (n-3/2)(n  +  1/2)  cos  (n  +  1/2)9  1+  (-i)n  d2n  r(n+1)  (A  -  3.d) 
(n-l)2Cos(n-l)0    -    (n+1)2  cos  (n+l)«l  \ 
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2_  NOOFM 

a  x     = 
aaar 


NOOFM      / 

£  |("1)(n~1)d2n-l  (n+l/2)  r(n_l/2)     |(n-3/2)sin(n-3/2)« 

n=l,.2,...(  L 


-  (n-3/2)  sin  (n+l/2)J   +    (-l)n  (n+l)  d2n  rn  (A  -  3.e) 

|(n-l)  sin  (n-l)#   -    (n+l)  sin  (n+l)* 

To  solve  the  problem,  it  was  assumed  that  the  stresses  along  the 
boundary  of  the  free  body  of  the  beam  can  be  closely  approximated  from  the 
conditions  of  the  uncracked  beam.      According  to  this  assumption,  we  have 
the  following  boundary  conditions: 


Along  AB: 

T 

xy 

0 

a 

X 

»V. 

- 

Along  BC: 

T            = 
xy 

r' 

xy 

/op 
iA(adf/a)  sin  (ae/a)  cos  (ay/a) 


*eB         ^  a[a*    +  ^(adf/a)] 


/  w 

f/  \  Pa       /      a  sin(ae/a)sin(a//a) 

1(x>      '       neB    J      z— o : •     da 

0  [a6    +«A(adf/a)] 

/oo 
sin(a</a)cos(a//a)        da 

0  a3   +    ^(«df/a) 


where 

6 


2      „2 


f  (x)     =       d3    [cd  +  xd  -  2xc  -  cz  -  xz  ] 


12    r  d 


S<x>     =    ~&   [X    "   T    +    c] 
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Along  CD: 


Txy      =     0 


ax       =     0 


Along  DE: 


xy 


=     0 


x 


=     -    P/2  6B 


In  the  above  expressions,  one  has: 

11/3 
Eblb 


a    = 


_Efb 


(1  -  Vf)    EbIh 
Efb 


1/3 


for  plane  stress 


for  plane  strain 


,/f(adf/a)     =    gdf/a    +   sinh(qdf/a)  cosh  (adf/a) 
sin  h     (adf/a) 


B  =  width  of  the  beam 

b  =  B/2 

c  =  crack  length 

d  =  depth  of  the  beam 

dr  =  depth  of  the  foundation 

E,  =  modulus  of  elasticity  of  the  beam 

Ef  =  modulus  of  elasticity  of  the  foundation 

L  =  moment  of  inertia  of  the  beam 
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/  =      effective  length  of  the  beam 

P         =      total  applied  load 
2  e        =      length  of  the  loaded  area 


of  symmetry 


"■«■  TM,»"  II    '«      I 


P/2€B 

»■§     u   «   »    >    » 


IUf 


P/2eB 


Crack 


Figure  A-l.        Collocation  Grid  and  the  System 
of  Loading 


Having  mathematically  described  the  boundary  conditions,    the 
system  of  simultaneous  algebraic  equations  for  the  collocation  procedure 
is  set  up  in  the  matrix  form  as  follows: 
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Boundary 

r 

"-P/26B 

D  -  E 

ax 

0 

C  -  D 

ax 

"'x 

A  -  B 

°y 

= 

°y 

B  -  C 

Txy 

0 

D  -  E 

rxy 

0 

C  -  D 

Txy 

xxy 

B  -  C 

Txy 

0 

A  -  B 

Expressing  the  above  system  of  equations  in  terms  of  the  Williams 
stress  function  coefficient,  dj  ,    yields: 


1,1 


'2,1 


CNOOFXY,  1 


'1,2 


'2,2 


Cl,2  (NOOFM) 


'2,2  (NOOFM) 


CNOOFXY,  2  (NOOFM) 


d- 


■ 


1 2  (NOOFM) 


*2 


(NOOFXY) 


(A  -4) 


where  C^  a  represent  the  coefficients  of  d^,    a^  represent  the  numerical  values 
of  the  boundary  conditions  to  a  particular  ith  boundary  station,  2 (NOOFM)  is 
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the  number  of  di  terms,  and  NOOFXY  is  the  total  number  of  equations  used 
in  the  program. 

The  program  will  solve  for  the  d-  ,     i  =  1,2,   .   .   .   .  2(NOOFM). 
Because  solution  convergence  could  be  a  problem,  the  least-squares  optimiza- 
tion technique  was  employed  in  solving  the  system  of  equations.    This  technique 
and  its  computer  listing  is  described  in  the  IBM  computer  manual  "The  Scien- 
tific Subroutine  Package"  —  Subroutines  LLSQ   and    DLLSQ,  pp.   160  -  164. 
This  method  of  solution  has  an  advantage  in  that  the  number  of  equations  can 
exceed  the  number  of  unknowns:    i.e., 

2(NOOFM)  <   NOOFSY  (  A  -  5) 

In  the  program,  NOOFXY  is  calculated  from: 

NOOFXY  =  2  [ NOOFXL  +  NOOFXU  +  NOOFXS  +  NOOFXB j  -  4 

(A-6) 
where  NOOFXL,  NOOFXU,  NOOFXS,  and  NOOFXB   are  the  number  of  colloca- 
tion stations  under  the  load  on  the  top  boundary,  outside  the  load  on  the  top 
boundary,  on  the  side  boundary,  and  on  the  bottom  boundary  respectively     In 
view  of  Equation  (A-6),  Equation  (  A  -  5)  becomes: 

NOOFM  +  2  <  NOOFXL  +  NOOFXU  +  NOOFXS  +  NOOFXB       ( A  -  7  ) 

This  equation  gives  the  criterion  for  the  allowable  size  of  the  system  of 
equations  (  A  -  4  ). 

Once  the  values  of  d^  have  been  solved,  the  stress-intensity  factor 
is  then  calculated  from: 

Kx     =     -  V^tt        di  (A  -  8) 
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and  the  corresponding  normalized  value  K'^  , 

K-, 


where 


K'l  = 

amax  Vd~ 

cr 

max 

6Pa2             1 

7T€Bd2            J 

0 

/•  00 

'              sin(ae/a) 
a      +tA(adf/a) 

da 

(A  -  9) 


(A-  10) 


Note  should  be  taken  of  the  fact  that  the  optimum  number  of  the  Williams 
stress  function  coefficients,  dj,  depends  somewhat  on  the  geometrical  and 
loading  conditions  of  the  beam.     Therefore,  a  sensitivity  analysis  similar 
to  the  one  shown  in  Figure  A-2  should  be  done  to  determine,  at  the  least,  the 
range  of  di  terms    (  =  2(NOOFM)  )   within  which  lies  the  maximum  value  of 
the  stress-intensity  factor.    Although  it  was  noted  that  the  optimum  value  of 
dj  terms  varies  slightly  with  the  crack  length,  only  one  set  of  sensitivity 
analysis  is  practically  sufficient  with  regard  to  the  economical  viewpoint; 
the  average  values  of  the  stress-intensity  factors  for  other  crack  lengths  of 
the  system  computed  from  three  or  more  values  of  d^  terms  within  the  opti- 
mum range  should  give  sufficiently  accurate  results  for  practical  purposes. 
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Example  of  Sensitivity  Analysis 
Input  data?: 


Beam:         Effective  length        =   3  inches 
Width   B  =1  inch 

Depth   d  =3  inches 

Young's  modulus  E^=    106  psi 


Foundation: 


=   semi-infinite 


Depth     df 

Young's  modulus  Ef  =  10**  psi 

Others:       Total  applied  load  P    =  10  lbs. 
Loaded  area's  length  =      1  inch 

Station  spacing  =  1/8  inch 

Crack  length  =  0.6  inch 


125 


100 


Number  of  6t  terms 


Figure  A-2.  Value  of  Normalized  Stress-Intensity  Factor 
Against  Number  of  Williams  Stress  Function 
Coefficients,  dj. 
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IIA.      USE  OF  THE   PROGRAM 

A.      Input  Data 

The  first  step  in  the  analysis  is  to  select  a  suitable  effective  length,  /  , 
of  the  beam  and  a  proper  station  spacing,  s.      It  is  recommended  that  the 
effective  length  be  greater  than  or  equal  to  the  beam  depth,  d,  and  the  sta- 
tion spacing  be  not  greater  than  1/8  of  the  beam  depth.    An  equal  spacing  be- 
tween each  boundary  station  is  assumed  in  the  program.    The  next  step  is  to 
choose  an  initial  value  of  NOOFM  (one-half  of  the  number  of  the  Williams 
stress  function  coefficients).    After  all  the  computations  corresponding  to 
this  value  of  NOOFM  and  a  specified  crack  length  are  completed,  the  value 
of  NOOFM  will  automatically  be  increased  by  an  increment  of  5  and  the  com- 
putations are  repeated  for  that  particular  crack  length  and  the  new  NOOFM. 
This  process  will  continue  through  the  last  NOOFM  specified.     The  value  of 
NOOFM  will  then  be  automatically  set  to  the  initial  value  again,  and  the  exe- 
cution continues  with  a  new  crack  length.     The  program  will  stop  executing 
when  all  the  computations  for  the  last  specified  crack  length  have  been  carried 
out. 

The  following  is  a  sequence  of  punched  cards  which  numerically  define 
all  the  parameters  needed  by  the  program. 

1.     Method  of  Analysis  for  Boundary  Conditions  (I  1) 
This  is  the  first  card  in  the  data  deck  to  give  the 

program  a  control,  whether  the  plane  stress  or  plane 
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strain  analysis  is  to  be  used  in  the  establishment  of  the 
boundary  conditions.    If  the  number  in  column  1  is 

0  the  plane  stress  analysis  is  to  be  used 

1  the  plane  strain  analysis  is  to  be  used. 

2.     Job  Title  (18A4) 

This  card  will  give  the  descriptive  identification  for 
the  job.    It  is  the  second  card  in  the  data  deck. 


3.     Geometrical  Properties,  Material  Properties,  and  Control 
Information  Cards 


The  third  card  through  the  seventh  card  will  give 
the  following  information: 

Card  Number  Format         Information  Columns 

3  3D24.16         Effective  length,  /  ,  of  the  beam  1-24 

Depth,  d,  of  the  beam  25-48 

Width,  B,  of  the  beam  49  -  72 

4  3D24.16         Young's  modulus,    Eb,  of  the  beam  1-24 

Young's  modulus,   Ef,  of  the  25-48 

foundation 

Depth  of  the  foundation,  df,  49  -  72 

(0  for  semi-infinite) 

5  3D24.16         Poisson's  ratio,  vf,  of  the  1-24 

foundation 

One-half  of  the  loaded  area's  25-48 

length,  € 

Total  applied  load,  P  49-72 

6  3D24.16        Initial  crack  length  1-24 

Crack  increment  required  25  -  48 

Final  crack  length  required  49  -  72 
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(2110,  L10,D24.16)        Initial  NOOFM  1-10 

Final  NOOFM  required  11-20 

Logical  constant,  PRISCK  21-30 

(.TRUE.,  or  .  FALSE.) 

Station  spacing,     s.  31-54 


B.      Output  Information 

The  following  information  is  developed  and  printed  by  the  program: 

1.  Reprint  of  geometrical  and  material  properties. 

2.  Stress-intensity  factor,  K-. . 

3.  Normalized  stress-intensity  factor,  K'^. 

4.  Maximum  tensile  stress  in  the  uncracked  beam, 
based  on  the  infinite  beam  on  elastic  solid,  crmr>v. 

5.  Fundamental  length  of  the  beam,    a. 

6.  Relative  crack  length,  c/d. 
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C * ***  THIS  PROGRAM  CALCULATES  THE  'Kl'  STRES S-INTENSI T Y_F ACTOR 

C****         ~  OE FINED  BY  : 

C****  Kl  =  -2.506628274631DO*SOL( 1) 

C**** WFTFlTE~5UrnT  TT""TRFnrTTT5n~TTJFFrn7JTT^^ "~ 

C#***  STRESS  FUNCTION 

C****  NCODE  rsrO"'RJFT"PDlTfE"'5Tft'E'55v~r^ FOR"~PLANE  STRAIfsT" 

C****  COND  IS  THE  BEAM'S  EFFECTIVE  LENGTH 

C****  CONH  IS  THE  BEAM'S  DEPTH 

C****  CONW  IS  THE  BEAM'S  WIDTH 

ZW*WZ  UME    IS    I  HE    BbAM»5~Y"OUNG''S   WOTJUDJT 

C****  CONEF  IS  THE  FOUNDATION'S  YOUNG'S  MODULUS 

C**~**  CONFD  IS  THE  FOUNDATION'S  DEPTH 

C****  CONVF  IS  THE  FOUNDATION'S  POISSON'S  RATIO 

C***#  CONA~  I~S~~ ONE-HALF    OF    THE    LOA  DE D   A R E A"» "S~L E NGT H" ~ 

C****  CONP  IS  THE  TOTAL  APPLIED  LOAD 

C*?*'*  CI  IS  THE  "CRACK  LENGTH  REQUIRED  W    THE  'COMPUTaT ION 

C****  CINC  IS  THE  INCREMENT  OF  THE  CRACK  LENGTH 

O***  C1END  IS  THE  FINAL  CRACK  LENGTH   REQUIRED 

C#***  MOOFM  IS  ONE-HALF  OF  THE  NUMBER  OF  WILLIAMS  STRESS  FUNCTION 

C****  COEFFICIENTS  REQUIRED  IN  THE  COMPUTATION 

C****  LNOOFM  IS  THE  MAXIMUM  VALUE  OF  " NOOFM  REQUIRED 

C****  PRISTJK  TS  THE'  LOGICAXCrjN$TANS"TANTr"nTRUEV»  INDICATES  THAT 

C****  AN  ACCURACY  CHECK  OF  THE  SOLUTION  WILL  BE  MADE 

C***#  SPACE  IS  THE  REQUIRED  SPACING  BETWEEN  EACH  BOUNDARY  STATION 

r*  <Ju  .»*.  »u  o* 
[^  ■*.*  'o  *r  'f 

IMPLICIT    REAL*"8    <A-H,0-Z) 

REAL*8  LAGUER 

REAL*4  NAM'ETIfiT 

REAL*4  FXL,  FXS,  FXB 

LOGICAL   PRISCK,ZEROCK,SOLUGK,BASEK 

COMMON  /DOU6/  PRIS»  B,   D  ,  CONE,  COR,  CORR,  COMA,  CONP,  CONH, 
*COND, C  1 ,  CONEF  ,  CONW,  CONFD,  C  1  END,  C  I  NC  ,  CONVF  ,  SPACE  ,  AM'AX  ,  S  I  GMAX  , 
*CONI,C0NSTA,DELXL,DELXU,DELXS,DFLXB 

COMMON  /SING/  NOOFM,  NOOFXU,  NOOFXS,  NOUFXB,  NOOFXL »  LNOOFM,  IARG, 

*  JARG,KARG,LARG,NAMF,NOOFI , NCODE , NOOFTO, MOOFXT , NOOFT 1 , 

*  NOOFTS,MOOFT2,NOOFTX,NOOFXX, NOOFY1 , NOOFYY, NOOFB1 , NOOFC 1 , NOOFOl , 

*  NOOFFl,NOOFGl,NOOFXY,NOOFS 
COMMON  /LOG/  PRISCK  ,  BASEK 

DIMENSION   X (001 00 ) ,Y( OOIOO) ,R( OOIOO) , PHI ( OOIOO) ,Vv( 00200)  , SOL (010 
*0),COFF(0200,.100),C0EFCK(200,100),CALC(200),CADD(2O0) , COM (200) , 
*CONCHK< 0200) , I P I  V ( 0100) 
EXTERNAL   FUNCTN 
READ(  5,  U6)NC0DE 
116  FORMAT(  I  1  ) 
3  CONTINUE 

BASEK  =. FALSE. 

READ( 5,8  01,EN0=1001 )  ( NAME (  I  ) ,  I  =  1 , 1 8 ) 
801  FORMAT ( 18*4) 
C 

C      READ  I'N  CONSTANTS 
C 

PPIS  =  1. 00-14 

RE AO ( 5 , 800 )COND, CONH , CONW, COME, CONEF, CONFD, CONVF , CUNA , CONP,C 1 , 
*C  INC, CI  END 
IF ( CONFD. GT.O. 1 )  BASEK  =  .TRUE. 
800  FORMAT ( 3024. 16) 

READ(5,802)  NOOFM,  LNOOFM,  PRISCK,  SPACE 
NOOFI   =   NOOFM 
*02  FORMAT (  21  10, L10, 024. 16) 

FXL  =  COMA/SPACE  +  1,05  191 


IMOO.FXL  =  IFIX(FXL) 
EXB  =  COND/SPACE  +  1.05 
NOOFXB  =  IFTX(FXB)' 
NOOFXU  =  MOOFXB  -  NOOFXL 
FXS  =  CONH/SPACE  +  1.05 
NOOFXS  =  IFIX(FXS) 
OELXL  =  COMA/(NOdFXL-r) 
DFLXU  =  (COND  -  CONA)/NOOFXO 
OFLXS  =  CONH/(NOOFXS  -~TJ 

OELXB  =  C0N0/(N_00^XB_-_1J_ 

M OOF TO  =  NOOFXL  +  1 

MOOFXT  =  MOOFXL  +  NOOFXO 

MOOFTi;  =  NOOFXT  +  ~T 

MOOFTS  =  NOOFXT  +J^00FXS_-_1 

N00FT2  =  MOOFTS  +  1 

NOOFTX  =  MOOFTS  +  MOOFXB  -  1 


M 0 0 F X X '"=" '"N "0 0 F XT  ~+~N QTJFXB- 

MOOFY1  =  MOOFXX  +  1 

NOOFYY  =  MOOFXX  +  NOOFXS"" 

MOOFfil  =  NOOFYY  +  1  

NOOFCl  "•=  MOOFYY  +  NOOFXT  -  1 

MOOF01  =  NOOFCl  +  1 

MOOFF1  =  MOOFCl  +  NOTJFXS  -2 

M00FG1  =  MOOFF1  +  1 

NDOFXY  =  MOOFYY  +  NOOFTX  -2 

MI  TWIT  =  MOOFXL  +  NOOFXO+  NOOFXS  +  NOOFXB 

IARG  =  MOO FT X 

LARG  =  MOOFXY 

COMB  =  COMW/2. ODO 

COM  I  =  C0MW*C0MH**3/12.0_D0 

CONSTA  =  (C0NE*C0Nl'7CdNEF7C0N~8)'**(  1.000/3.000) 

IFTMCOOE.  EO.  1  )COMSTA=(  1 . 0D0-C0NVF**2  )**'(  1 .  000/ 3.000  )_*CJJNSTA_ 
C 

C***4     CHECK  FOR  SEM I- I NF I NI TE  DEPTH  FOUNDATION. 
C 

\     IF(BASFK)  D  =  CONFD/CONSTA 

D  =  10000. ODO 

B  =  CONA/CONSTA 
:  AMAX  =  LAGUERI  _",  B  ,  0.  ODO  , IT,  FUNcTn")  " 

SIGMAX  =   6.0D0*AMAX/(C0NH-v*2*C0NW) 
9000  CONTINUE 

IF ( NOOFM. GT.  LNOOFN)  GO  TO  9999 

IF  (NOOFM. EO.'NOOFl")  GO  "TO"  9997 

GO  TO  999B 
9999  CI  =  CI  '  +  "C'INC 

NOOFM  =  MOOFI 

C1EMD  =  C1EMD+0.1D-3 

IF(C1.GT.  CI  END)  GO  TO  3 
9997  CALL  COORD  (X,Y,R,PHI) 
999B  CONTINUE 

KARG  =  NOOFM* 2 

JARG  =  2*KARG 

NOOFS  =  KARG 

CALL  STRESK        (  X,  Y,  R,  PH  I ,  V9,  SOL  ,  COEF,  COEFCK,CALC  ,-CADO,  CON, 
*  CONCHK, IP  IV ) 

NOOFM  =  NOOFM  +  5 

GO  TO  9000 
1001  STOP 

END 

SUBROUTINE  COORD   (X,Y,R,PHI) 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*4  NAME( 18 ) 
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COMMON  /DOUB/  PRISt  B,   0  ,  COMB,  COR,  CORK,  COMA,  CONP,  CONH, 
*COND, Cl,CONEF, CONW, CONFD, C 1END, C I NC , CONVF , SPACE, AM AX , SI GMAX , 
*C0NI,C0NSTA,DELXL,DELXU,DELXS,DELXB 

COMMON)  /SING/  NOOFM,  NOOFXU,  NOOFXS,  NOOFXB,  NOOFXL,  LNOOFM,  IARG, 

*  J ARG, KARG, L ARG, NAME, NOOF I , NCODE, NOOFTO, NOOFXT , NOOFTl , 

*  NOOFTS  , NOOFT2,  NOOFTX,  NOUFXX ,  MOOF  Y  1 ,  NOOFYY,  NOOFBl  ,  NOOFC  1 ,  NOOFD  1 , 

*  MOOF F 1 , NOOFG 1 , MOOF X  Y , NO( I F S 

DIMENSION  X(IARG),  YdARG),  R(IARG),  PH'KIARG) 
C 

C         GENERATION  OF  X  S  V 
C 

Y(  1 )  =  0. 

X(  1  )  =  CONH  -  CI 

DO  803  1=2, NOOFXL 

Y( I )  =  Y(I-  1)  +  DELXL 

803  X( I )  =  X ( 1-1 ) 
Y(NOOFX|_  )  =  CONA 

DO  863  I=NOOFTO,NOOFXT 
Y(  I  )  =  Y(  I-  1)  +  DELXO 
863  X( I )  =  X( 1-1 ) 

X(NOOFXT )  =  CONH  -  CI 
Y(NOOFXT)  =  COND 
DO  804  I=N00FT1,N00FTS 
Y( I )  =  Y( 1-1) 

804  X(  I  )  =  X( I-  1 )  -  DELXS 
Y(NOOFTS)  =  COND 
X(NOOFTS)  =  -  CI 

DO  805  I=M00FT2, NOOFTX 
Yd)  =  Yd-  1)  -  DELXB 

805  X( I )  =  X( 1-1 ) 
Y( NOOFTX)  =  0. 

C 

C         CALCULATION  OF  POLOR   COORDINATES 

C, 

DO  "5  1  =  1,  NOOFTX 

R(  I )=DSORT( X(  I  )**2  +  Y(  I  )**2) 
5  PHId  )=DATAN2(Y(  I  ),X(  I  )  ) 
C 

RETURN 
C 

END' ~-" 

C 

SUBROUTINE  STRESK  ( X , Y , R , PH I , V9 , SOL ,COEF , COEFCK , C ALC ,C ADD, CON , 

*  CONCHK, IPIV) 

TMPCIC  I T  "REAl^B    (  A-H ,  0-T  J 
REAL*8    LAGUER 

"  "READ- 4"        STAKT3G  ) ,  NTTWETTHl " ~ — 

LOGICAL       PRISCK,ZEROCK,SOLUCK,BASEK 
"COMMON    /DOUB/    PR!  ST  "B;    "D    t^COME,    COR,    CORR",    CONA",    CONP,    CONH, 
*COND, C 1 , CONEF , CONW, CONFD, C 1  END, C I NC , CONVF , SPACE, AMAX , S I GMAX  , 
".  *CON I ,  C ON S1XTDE L XT~,T)E L X"U ,  D'ElTsTDE  L  X  B 

COMMON  /SING/  NOOFM,  NOOFXU,  NOOFXS,  NOOFXB,  NOOFXL,  LNOOFM,  IARG, 

""  ~~*-jmG7Km&;  i  a  r  g',  n  ame  ,  n  o  u  ft  ^codttwoftct,  no  o  f  xttno  tt  fit; 

*  NOOFTS,NOOFT2,NOOFTX,NOOFXX,NOOFY1,NOOFYY,NOOFB1,NOOFC1,NOO FD U 

*  N00FF1,N00FG1,N00FXY,N00FS 
COMMON  /LOG/  PRISCK  ,  BASEK 


DIMENSION  X(IARG),  Y(IARG),  R(IARG),  PHIdARG),  V9UARG),  SOL(KARG 

2  ) ,  CO  E  F  (  L  ARG,  KARG  ) ,  COEFCK  (  LARG,  KARG  )  ,  CALX  (  LARGJ.»C_A_DD_Lk  ARG)  , 

TT"ON  (  L ARG  ) ',  CONCHK (  L ARG  )  ,  I  PI  V  (  KARG  ) ,  F0UNDT4  ) 

DATA  STAR/30* '****'/, FOUND ( 1)/ 'SEMI »/,F0UND(2)/ '- INF ' /, FOUND ( 3 ) / 
*  •  INIT'/,F0UND(4)/ »E'/ 

EXTERNAL    FUNCTN  193 


ZEROCK  =  .TRUE. 

SOLUCK  =  .TRUE. 

IF(Cl.EO.O.O)   ZEROCK  =  .FALSE. 

IF(Cl.EO.O.O)   PRISCK  =  .FALSE. 

IF(  (N00FM  +  2)  .GE.  (  NOOFXL+NOOFXU+NOOF.XS+NOOFXB )  ) 

IF( ( NOOFM+2 ).GE. ( NOOFXL+MOOFXU+NOOFXS+NOOFX8 ) ) 

IF( ( NOOFM+2 ).GE. ( NOOFXL+NOOFXU+NOOFXS+NOOFXB ) ) 

IF(Cl.EO.O.O)   GO  TO  8000 


ZEROING  BLANK  COEFFICIENTS 


PRISCK=. FALSE. 
SOLUCK*. FALSE. 
GO  TO  8000 


C 
C 

c 


DO  22  I  =  ltNOOFXY 
CONCHM  I  )  =  0. 
CON( I )  =0. 
DO  22  J=1,N00FS 
COEFCKI I, J)  =  0. 
22  COEF( I, J)  =  0. 

CALCULATION  OF  SIGMA  XX  COEFFICIENTS 


DO  15 
N  =  I 
IF(  I. 
IF(  I. 
IF(  I. 
0  =  1.  D 
DO  15 
IF(  J- 
C0R  =  1 
GO  TO 
COR  =  - 
CORR  = 
ESIN  = 
ECOS  = 
IF(DA 
IF  (DA 
CALL 
COEF( 
**ECOS 
CALL 
COEF( 
**EC'OS 
15  0=0+1 


I=1,NOOFXX 


\ 


N  =  I  +  NOUFX'S  -  2 

CON(I)  =  -CONP/( 2.0D0*C0NA*C0NW) 

CON(  I  )  =-LAGUER( 3, B , Y ( N ) /CONST A , D , FUMCTM ) /CONW 


80 

81 
201 


GT.MOOFXT) 
LE.NOOFXL  ) 
GT.MOOFXT) 
0 

J=1,NQ0FM 
(  J/ 2* 2")  )     80T80»8l 
.DO 

201 
l.DO 
-COR 

DSIN(  PHUN)  ) 

DCOS( 'PHI(N')  ) 
BS(ESIN).LT.PRIS)       ESIN=0.     , 
BS(ECOSr.Lf.PRIS)       ECOS=0. 
STR  E  S  1_   (0,PH  I  (  N  UO 1,D2  L0  3  ?  D4  ,_D5_,  R  (N  )  ) 
I,2*j-1)    =         D1*ESIN*E;SIN      +D2*EC0S*EC0S/R  ( 
/R(N)       +D5*ECOS*ECQS/R(N)**2      -      D4*2.D0*EC 
STRES2 '  ("Q,  PH  I  (  N  )  T  Dl ,  D2,  D3,  D4,  D5~,~R( "N  )  ) 
I,2*J)     =  D1*ESIN*ESIN       +D2*EC0S*EC0S/R ( 


N) 
0S: 


+    D3*2 
■ESIM/R( 


.[)0*ESIN 
N)**2 


/R(N) 

.DO 


N) 

+  D 5  *  E COS*ECOS/R( IM )  *  *  2      -  "    04*  2 . D 0*  E COS  * 


+  03*2.D0*ESIN 
ESIN/R(N)**2 


C 
C 
C 


CALCULATION  OF  SIGMA  YY  COEFFICIENTS 


C  =  _Y(N 0 0 F X T  ) /CONST A_ 

CONM  =  LAGUER( 1 , 8, C, D, FUNCTN j 

_D0  17  I=N00FY1,N00FYY 

~N 


=     I     -    NOOFXB    -    1 
CON( I )=-CONM#(X(N)-CONH/2.D0+Cl )/CONI 


84 


85 
203 


0=1. DO" 

DO    17.  J=1,N00FM 

I F  (  J  -  (  J  77*271    84,8  4,  8  5" 

C0R=1.D0 

GO"  TO    203 

>C0R  =  -1.D0 

CORR  =  -C'OR  """ 

ESIN=       DSIN(PHKN)) 


~ECOS  =       DCO~S(  PHI(N)  ) 
IF(DABS( ESIN).LT.PRIS) 


ESIN=0. 


% 


T9T 


TF  ( DA  B  S  T  E  CO  5T.TT.WrST     EXUS  ^OT 

CALL    STRES1    (  Q,  PH  I  (  N  )  ,  Dl ,  D2,  03,  D4,  D5  ,  R  (  N  j  ) 


C0EF(  If2*J-lT    =      "  01*EC"0S*EC0S       +D2*ES  IN*ES  I'N/RTN  )"      -2.D0*D3*ESTn* 

*EC0S/R(N)       +    2.D0*D4*EC0S*ESIN/R('N)**2       +      D5*ES IN*ES IN/R ( N ) **2 

CALX  "STRTS2~TD\~F^rriTr,XTT0T,Tn-7D^DTrRTNTl ~       " 

C0EF(I,2*J)  =      D1*E_C0S*EC0S   +D2*ES IN*ES IN/R ( N)   "2. D0*D3»ESlN* 

*  EC  OS  /  R  (  N  ) "  "  +  2  .  DO*  04*  EC  OS*  ES I N  ?R  (  N  )  **2   + "  ~"D5*  E  S I  N*E  S  IN/  R  (  N ) "**  2 
17  0=0+1.00 

C 

C      CALCULATION  OF  SIGMA  XY  COEFFICIENTS 

C  — -  - 

CONG  =  -LAGUER( 2, B, C, D, FUNCTN ) /CONW 

DO  16  I=NOOFB1,NOOF~XY 

N  =  I  -  NOOFYY  +1 

IF(M.LE.NOOFXT"VORT  N.GE.NOOFTS)  GO  TO  550 

C0N( I )  =  6*0DO*C0NQ/CONH**3*(Cl*C0NH+X(N)*C0NH-2.0D0*X(N)*Cl-Cl**2 

*  -X(N)**2) 
550  0  =  1.00 

DO  16  J=1,N00FM 

IF( J-( J/2*2) )  82,82,83 

82  C0R=1.D0 
GO  TO  202 

83  C0R=-1.D0 
202  CORR=-COR 


FC'OS  =   DCD  ST2 .  D '0* P Hi  ( W)  ) 
ESIN  =   DSIN(PHKN)  ) 

ECOS  =   DCOS(PHKN)) 

IF ( DABS ( ESIN).LT. PRIS)   ESIN=0. 

IF(DABS!FCOS).LT.PRISJ  ""  FC0S=0. 

IF(DARS(ECOS).LT.PRIS)   ECOS=0. 

CALL  STRES1  ( 0 , PH I( N ) , 01 , 02 , 03, 04r D5 , R ( N ) ) 

C0EF( I,2*J-1)  =   -D1*EC0S*ESIN    -D3*FC0S/R ( N  )   +p5*EC0S*ESIN/R ( N) 
***2    +D2*EC0S*ESIN/R(  N)   +D4*FC0S/R(  N  )**2 

CALL  STRES2  ( 0, PH I ( N ) , 01 , D2 , 03 , 04, 05 , R ( N ) ) 

C0EF!T,2*JJ  =     =D1*XCCTS*E"SIN    -D3*FC0S/R  ( N )   +D5*BC0S*ESIN7R'(N) 
***2    +D2*EC0S*ES  I N/R  (  N  )   +04=:=  FCOS/R  (  N  )  **2 
16  0=0+1.00 
C 

C****     SOLVING  FOR  SOLUTIONS 
C 

00  130  I  =  1,N00FXY 

COMCHK( I )  =  CON( I ) 

DO  130  J  =  1,N00FS 
130  COEFCK! If J)  =  COEF( I, J) 

CALL  DLLSQ(COEFCK,CONCHK,NOOFXY,NOOESt 1, SOL , IPIV, 1 .E-72t IER,V9) 

SFFK1  =  -2.50662827463100*SOL( 1) 

YY  =  SEFK1/(SIGMAX*DSQRT(C0NH) ) 

CH  =  Cl/CONH 
C 

C         OUTPUT 
C 


8000  WRITE(6t810)    ( STAR ( I ) , I =1 , 30 ) 
810  FORMAT! IH1/1H1  //   30A4) 

WRITE! 6, 873)  ( NAME (  I  ) ,  I  =  1 ,  1 8 ) 

873    FORMAT!  //30X,  'STRESS       INTENSITY       FACTOR       A    N 

*    A    L    Y    S    I    S»///25X,  'CASE    :        ',18A4) 
IF<NCOOE    .LE.    0)WRITE(6,114) 
IF     (iMCODE    .  EO.     1)     WRITE!  6,  115)C0NVF 

114  FORMAT! //30X, 'PLANE  STRESS  SOLUTION.') 

115  FORMAT! //30X,  'PLANE  STRAIN  SOLI.JT  I  ON.  • ,  21X  ,  •  FOUNDAT  I  ON  •  '  S  POISSON1' 
*S  RATIO  =' ,F7.4) 

C0N2A  =  2.0D0*C0f\!A  195   o- 


8  72 


6719 


671* 


*76 


IF( «AS 
WRITE ( 

*co-vp,c 

30  TO 
31?  WRITE ( 

*CPMPf.C 
313  COM TIN 
871  FORMAT 

*F10. 1, 

*  'CRACK 

*  :  •  /  /  3 
-"•OTHER 
*/30X,  » 

FORMAT 

*  F 1 0  .  1 , 

*  'CRACK 

*  :  '  /  /  3 
* 'OTHER 
*/30X,  ' 

IFTSOL 

WRITE ( 

FORMAT 

EXECUT 

IFIZER 

W  R  I T  E  ( 

FORMAT 

♦LENGTH 

GO  TO 

879  IF(  IE.R 

WRITE( 

90  FORMAT 

* L  Y  ILL 

8  78  WRITE( 

299  FORMAT 

'■   WRITE( 

8  13'  FORMAT 

fe//30X, 

3  F12.6 

WRITE'( 

IF( PRI 

GO  TO 

W  R  I  T  E  ( 

FORMAT 

DO  3  2 

C  A  L  C  (  I 

DO  3  2 

CALC(  I 

DO  526 

CADD(  I 

WRITE( 

♦NOOFBl 

296  FORMAT 
S///43X 

*  ,  1 4 ,  ' 

*  'FOOAT 
*LOCATI 

*  D  A  R  Y  •  , 
#/30X, • 

*  'FOOAT 

*  'FOOAT 
WRITF( 

297  FORMAT 


OF  BOUNDARY  STATIONS 


SPECIFIED  CRACK 


877 


293 
8107 


32 


5268 


EK)  GO  TO  312 

6.871  )  CONE,  CONH,COND.C1,CONW,  CONE  F,   ( FOUND (  I  )  ,  1  =  1 ,4)  ,    

0N2A,NOOF~S, SPACE 

313 

6. 872  )00NE,CONHtCOND,Cl,C0NW,CONEF,CONFD, 

0M2A,N00FS, SPACE  

OF 

(7/56X,  •  PARAMETERS  »//25X,  'BEAM   :  V/30X,  '  YOUNG ^  «S  MODULUS'  , 

18X,  'BEAM  ObPTH « , 5X, F10.2, /73X, 'BEAM  LENGTH'  ,4X , FLO. 2» /30X , 

LENGTH  • ,  IX,  Fib. 4-,  20X,  •  BEAM  WIDTH  »  ,  5X  ,  F  1  0.  2//25X  ,  '  FOUNDATION 
OX.,  '  YOUNG  •  '  S  MODULUS  '  ,  F 1  0.  1,1 8X  ,  '  DEPT H  '",10X",  4A4 ""  ,7/25X  , 
S   :  '//30X,  'LOAD',  9X.F10,.  1,20X,  'LOADED  AREA'  'S  LENGTH'  ,F8;4, 
NUMBER  OF  SOLU-T  I  ONS  '  ,  T4,"20X,  •  ST  AT  I  ON  SPAC"  I  NG  '  ,4X  ,  F8 .3 j 
( //56X,  'PARAMETERS '//25X,  'BEAM   :  '//30X,  'YOUNG' k S  MODULUS' , 
18X,  'REAM  DEPTH' ,5X, F 10.2,/ 73X,  'BEAM  LENGTH'  ,  4X  ,  F  10.  2  ,/30X,  " 

LENGTH',  IX,  F10.4,20X,  'BEAM  W  I  DTH  •  ,  5X  ,  F  1  0.  2//25X  , _'  FOUND  AT  I  0N_ 
OX , • YOUNG ' • S  MODULUS • , F10. 1 , 1 8X , 'DEPTH ' VlOX ,  F12 .5 , 7/2  5  X , 
S   :  '//30X,  'LOAD', 9X, F10. 1,20X,  'LOADED  AREA''S  LENGTH ', F8. 4, 
NUMBER  OF  SOLU.Tl  ONS  '  ,  I  4,  20X,  '  STAT  ION  SPAC  ING  •  ,  4X",  F8  .  3  ) 
OCK)   GO  TO  6718 
6,6719) 

(//15X,  •*  *  *  *  *  NOOFM  >  NUMBER 
ION  WAS  TERMINATED  *  *  *  *  *') 
OCK)   GO  TO  8  79 
6,876) 
(//   15X,  •*  *  *,  *  *  THEORY  IS  NOT  VALID  FOR 

—   EXECUTION  WAS  TERMINATED  *  *  *  *  *•) 
87  7 

. EO.OJGO  TO  878 
6,90) 

(///01X»*  *  *  *  *  SOLUTION  PROCEDURE  INDICATES  MATRIX 
-CONDITIONED  OR  RANK  LESS  THAN  NUMBER  OF  COLUMNS  *  *  * 
6,299)   SEFK1 

(   ///   "30X,  'STRESS  INTENSITY  FACTER  Kl  =',024.10) 
6,813)  SIGMAX,  YY,  CH,CONSTA 
(   /, 30X,  'MAXIMUM  STRESS  IN  UNCRACKED  BEAM  =  «,F12.6, 

•  NORMAL  I  ZED  K  =  ', F 1 2 . 6, //, 30X ,' CRACK/DEPTH  RATIO  =', 
, //30X, 'FUNDAMENTAL  LENGTH  OF  BEAM  =  «,F12.6) 
6,811)    ( STAR(  I  )  ,  1  =  1,30) 
SCK)  GO  TO  293 
1001 
6,8107  ) 
(  1H1// 

1=1,  NOOFXY 
)=0. 

J=l,  MOOFS 
)=  CALC(I)  +  SOL( J)*COEF( I, J) 
8  1=1,  NOOFXY 
)  =  COM( I )-CALC( I ) 

6,  29  6)  NOOFXL*NOOFTO,NOOFXT,NOOFT1,  NOOFXX,  N00FY1.,  NOOFYY, 
, MOOFC 1 , N00FD1 , NOOFF 1 , N00FG1, NOOFXY 

(///35X,'*  *  *  *  *  OOTPOT  CHECK  ON  SOLUTION  FIT  *  *  *  *  *  • 
,  'SUMMARY  OF  BOUNDARY  REQUIREMENTS • //30X , ' EQUAT I ONS     1  TO' 
SIGMA  XX  ON  TOP  BOUNDARY  —  LOCATIONS  UNDER  L0ADV30X, 
IONS  ',14,   '  TO', 14, '  SIGMA  XX  ON  TOP  BOUNDARY  —  REMAINING 
ONS ', /30X, /EQUATIONS  ',14,'  TO', 14,'  SIGMA  XX  ON  BOTTOM  BOON 
/30X, 'EQUATIONS  ',14,'  TO', 14,'  SIGMA  YY  ON  SIDE  BOUNDARY', 
EQUATIONS  ',14,'  TO', 14,'  SIGMA  XY  ON  TOP  BOUNDARY ',/30X , 


SEVERE 

:  *   *  •  ) 


( STAR( I ), 1=1,30) 
30A4) 


IONS  ',14,  •  TO' , 14,  'SIGMA  XY  ON 
TO' ,14,  'S  IGMA  XY  ON 


IONS  ', 14, • 

6,297) 

(  //  14X,  'EQUAT 10 


SIDE  BOUNDARY' ,/30X, 
BOTTOM  BOUNDARY' ) 


NO.  '  ,  13X,  'THEORETICAL  VALUE'  ,  12X, 'CALCUAT 
1 96  ' 


*ED  VALUE •  ,  15X,  "DIFFERENCE »   ) 
WRITE( 6,291  )   (  I,CON(  I  ),CALC(  I  ),CADD(  I  )  ,  I  =  1 ,  NGQFXY  ) 
291  FORMAT (  17X,  13, 10X, 302  8. 16) 

WRITE! 6, 8  11 )    ( STAR ( I ), 1=1,30) 
811  FORMAT*///    30A4) 
1001  CONTINUE 

RETURN 
END 


SU6R0U 
IMPLIC 
COMMON 

*COND,C 

*CONI,C 
ECOS  = 
FCOS  = 
ESIM  = 
FSIM  = 
IF(OAB 
IFIDAB 
IF (DAB 
IF (DAB 
02PDR2 

«RR«( 
DPDR=D 
D2PDRP 

*(+1.00 
DPDP=D 
D2PDP2 

*D0)*( + 


TIME  S 

IT  REA 

/DOUB 

1,C0!ME 

ONSTA, 

DCOS( 

OCOS( 

DSIN( 

DSIM( 

S  (  E  S  I  N 

S( ECOS 

S( FCOS 

S( FSIM 

=  (  -e 

S-1.5D 
2PDR2* 
=( (S-l 
)**(S- 
2PDRP* 
=  (  (  S 
l.DO)* 


TR 

L* 

/ 

F, 

DE 

(  S 

(S 

(  S 

(S 

). 

). 

). 

). 

CO 

0) 

RR 

.5 

1. 

RR 

-1 

*( 


ESI 

8  (  A 

PRIS 

CONW 

LXL, 

-1.5 

+  0.5 

-1.5 

+  0.5 

LT.P 

LT.P 

LT.P 

LT.P 

S  +( 

*  (  + 1 

/(  S- 

00)* 

DO)* 

/(  S  + 

.  500 

S-l. 


(S,PH 

-h,0- 

,  B, 

,CONF 

DELXU 

DO)*P 

DO)*P 

00)*P 

DO)*P 

RIS) 

RIS) 

RIS) 

RIS) 

(2.  DO 

.  DO )  * 

.500) 

ESIN 

CORR 

.5D0) 

)**2* 

DO)*C 


EE,L)2PDR2,DPDR,D2PDRP,DPDP,D2PDP2,RR) 
Z) 

0  ,  CONE,  COR,  CORR,  COMA,  CONP,  CONH, 
D, CI  END, C INC , CONVF, SPACE , AMAX, S I GMAX, 
,DELXS,DELXB 
HEE  ) 
HEE) 
HEF) 
HFE) 

ESIN=0. 

ECOS=0. 

FCOS=0. 

FSIN=0. 
*S-3.D0)/( 2.D0*S+1.D0)  )*FCOS) 
*( S-1.DO)*CORR 


(S**2-.25D0) 


-( S- 1 . 5D0 ) *  FS I N ) * ( S+ . 5D0 ) *RR** ( S- . 5D0 ) * 


ECOS 
ORR 


- ( S-l . 500 1* ( S+. 500 ) *FCOS ) *RR** ( S+. 5 


RETORM 
END 


SUBROU 
IMPLIC 
COMMON 
=  COND,C 
=CONI,C 
FSIN  = 
ECOS  = 
ESIM  = 
FCOS  = 
IF(DAB 
IF(DAB 
IF(DAB 
IF(DAB 
D2PDR2 
DPDR=D 
D2PDRP 
**COR 
DPDP~=D 
D2PDP2 


TIME  STRES2  ( S, PHEE , 02PDR2 , OPDR , D2PDRP , DPDP , 02 POP  2 , RR ) 
IT  REAL*8  (A-H,0-Z) 

/DOOB/  PRIS,  B,   D  ,  CONE,  COR,  CORR,  COMA,  CONP,  CONH, 
1, CONEF, CONW, CONED, CI  END, C INC, CONVF, SPACE, AMAX, SI GMAX, 
ONSTA, DELXL, DELXU, DEL XS, DEL XB 

DSIN"( (S+1.D0)*PHEE) 

DCOS( ( S-1.D0)*PHEE) 

DSIN( ( S-1.D0)*PHEE) 

DCOS( ( S+1.D0)*PHEE) 
S(ESIN).LT.PRrS)   ~  ESIN  =  0. 
S(ECOS).LT.PRIS)       ECOS=0. 
STF5TN1  .  LT.PRTS T   "TS  TN=0  . 
S(FCOS).LT.PRIS)       FCOS=0. 

(-ECOS  +FC0S)*(S*S+S)*RR**(S-1.00)*(+1.00)* 
2PDR2*RR/S 

(  ( S-l .  DOT* ES I N  -( S  +  l . DO ) * FS I N ) * ( S  +  l . DO ) *RR* 


*S*COR 
*S*(+1.D0)**S 


i2  pdr  p^RRrrs+"iTDcn 

=   ( ( S-l. DO)** 2* ECOS 
*T**S*COR  " 


- ( S  +  l . DO  )**2*FC0S )*RR** ( S+ 1 . DO ) * (  +  1 . DO 


C 
C 


RETURN 

END  __   __ 

REAL  FUNCTION  FUNCTN*8 ( NFUNCT, B, C , D, X ) 
IMPLIC  IT  RE~AL*8  (  A-H  ,  O-Z  ) 

197 


S  =  D*X 

IF  ( S.  GT. HO. )  GO  TO  2 

PS  I  =  (S+OSINH(S)*DCOSH(S) )/(DSINH( S) 

GO  TO  3 

2  ps I  =  l.noo 

3  IF(MFUfMCT.EO.  1  )  FUNCTN= 
IF ( NFUNCT. EO. 2)  FUNCTKl  = 
IF(MFOMCT.EO. 3)  FUNCTN= 

*{'X*'(X#*3  +  PSI  )  ) 


0SINH( S) ) 


( OFXP( X)*OSIN( B*X )*DCOS(C*X ) ) /( X**3+PSI ) 
(  X*DEXP(  X)*DSIN(  B*X')*DS  IN(  C*X )  )/(  X**3+PSI  ) 
(  DEXP(  XH=P-SI*DSIN(  B*X)*DCOS(C*X)  )  / 


RETURN! 

END 

REAL  FUNCTION  LAGUER*8<  NFUNCT,  BB,  C,DD,  FUNCTN  ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  /DOOB/  PRIS,  B,   0  ,  COME,  COR,  CORR,  COMA,  CONP,  CONH, 
'COND'f  C  1 .  ,CONEF,  CONW  ,  CONFD,C  TEND,  C  I  NC, CONVF, SPACE ,  AM  AX  ,  S  IGMAX  , 
'CONI,CONSTAfDE;LXL,DELXU,DELXSrDELXB 


C 

C, 

C 


DIMENSION  1(32) j W(32) 

PRESET  Z  AND  W  ARRAYS, 


DATA  PI 
DATA  1/ 
1  ,.1072 
34922 
73581 
12763 
19855 
28862 
40145 
54333 
72687 
98829 


/3.14 
.4448 
44875 
13273 
26733 
69798 
86094 
10181 
71977 
72133 
62809 
54286 


159265 
936583 
381781 
021994 
186241 
674272 
033605 
632347 
153944 
339690 
&66271 
828397 


3589 
3267 
76D1 
5  D 1 , 
D 1 ,  . 
5D2, 
5D2, 
5  02, 
2D2, 
702, 
D2,. 
D2t. 


793D0 
0180- 
,.172 
.4616 
89829 
.1493 
.2263 
.3234 
.4450 
.5989 
80187 
11175 


1,. 2 345 26 109 5 196 18 54, .5  76  884629  30188643 
2A08  7764446454D1, .252  8336706425794901, 
4^6769749767401 , .5903958504174243901, 
40  924212  59601, . 1078 30 186 32 53997 202, 
113975552255702, .1729245433671531502,  " 
088901319677402, . 2 562 86 3 60 22459248 D2, 
662915396473702, .3610049480575197402, 
920799 57 5493802,. 49 2243V4987 30863902, 
2509 162 134018D2, . 65975377287935053D2  , '" 
4469779 13 52 D2, .88  7  3  5340417  8924002, 
13980979377003/ 


DAT 

1  .  1 

2  .3 

3  ,. 


4 
5 
6 
7 
8 
9 
I 
J 
K 
L 


,  . 
.6 
.2 
.3 
.  1 
.2 
.5 
.  1 
.2 
.4 


A  W/.109 
95903335 
176  09125 
98080330 
59345416 
35060222 
30589949 
23780165 
54213383 
05442967 
66129413 
91337549 
67151121 
51053619 


2183419523 
97288104, . 
091750700- 
66149551D- 
128686329D 
66258067D- 
1 89 13 36  ID - 
77292665D- 
339  38  234D- 
37880454D- 
03973594D- 
44542243D- 
9  24013700- 
38989742D- 


849 
129 
It. 
3,. 
-5, 
8,. 
10, 
13, 

167 
20, 
25, 
30, 
37, 
47/ 


7, 
98 
11 
21 
.7 
42 
.9 


.2104431079  3881323, .2  35  21322966984801, 
378628607176, .705786238657174420-1, 


918214834838557D-1, . 3738  81629461 1 5248D-2 
48649 18801 364 19D-3,. 39 20341 967987947 2D-4 
41 6404578667552D-6,. 7604567879 1207810-7, 

8138297104092890-9,  

799379  2  8  8 72  7094D- 12, 

17 1 8234434207 19D-1 5, _ 

1 1979  2  2  9  01 6  3618  6  D  - 1 8, 


2 
1 
1 

13386169421062  5630-41, 


.346982  5  866  37395  20-22, 
41856054546303690-2  7, 
192  248760098  2  2240-33, 


SUM    =0.0 
DO    5    J  -■=    1,32 
5    SOM    =    SUM+W( J )*FUNCTN(NFUNCT, BB,C,DD,Z( J ) ) 


AREA    =    2.0D0*C0IMA 

I F ( NFUNCT. EQ. 1 )  LAGUER" 

I.F  (  NFUNCT.  EO.  2)  LAGOER 
IF(MF0NCT.E0.3)  LAGOER 


2.0D0*C0NP*C0NSTA**2*SUM/PI/AREA 
-2.0D0*C0NP*C0NSTA*SUM/P  I_/AREA___ 

2V0D0*C0NP*SUM/PT/AREA '" 


RETURN 
END 


198 


SUBROUTINE    DLLSO(  A»B_,  M,  N,  L,  X,  IP_IV,EPS,  IER,AUX) 

IMPLICIT  "REAL*  8    ( A-H,0-Z )  " 

REAL*4    EPS 

D I M  ENS  TO  WAT  1  J,B(1),X(  1  nTPTVT  I  F,  AUX  (  IT" 

1  PIV=O.DO 
IEND=0 

DO    4    K=ltN 

IPIV{K)=K 

H=0.00 

IST=IEND+T 

IEMD=IEND+M 

DO    2     I=ISTrIEND 

2  H  =  H  +  A(  I  )*A<  I  ) 
AUXfK  )=H 
IF(H-PIV)4,4,3 

3  PIV=H 
KPIV=K 

4  CONTINUE 

5  SIG=DSORT(PIV) 
TOL=SIG*ABS( EPS) 
LM=L*M 

IST=-M 

DO    21    K=1,N 

!ST=IST+M+1 

IEND=IST+M-K 

I=KPIV-K 

IF( I )8,H, 6 

6  H=AUX(K) 

AUX( K )=AUX< KPIV) 
AUX(KPIV)=H 
ID=I*M 

DO  7  1  =  1  ST, I  END 
J=I+ID 
H  =  A(  I) 
A(I }=A( J) 
A( J)=H 

IF(K-1  )  1.1,  11,9 
SIG=O.DO 

DO  10  1  =  1  ST, I  END 
SIG=SIG+A(  I )*A(  I  ) 
SIG=DSORT( SIG) 
IF(SIG-TOL  )  3  2, 32,  11 
H=A( 1ST) 
IF(H) 12,  13,13 
SIG=-SIG 

IPIV(KPIV)=IPIV(K) 
IPIV( K)=KPIV 
BETA=H+SIG 
A( IST)=BETA 
BETA=l.DO/( SIG*BETA) 
J=M+K 

AIJX  (  J  )=-SIG 
IF(K-M) 14, 19, 19 
14  PIV-O.DO 
10  =  0 
JST=K+1 
KPIV=JST 
DO  18  J=JST,N 
ID=ID+M 
H=0.00 
DO  15  I=IST,IEND  199 


DLLS  630 


7 
8 
9 

10 


11 

12 
13 


DLLS  730 

DLLS  740 

DLLS  750 

DLLS  760 

DLLS  770 

DLLS  780 

DLLS  790 

DLLS  800 

DLLS  810 

DLLS  820 

DLLS  830 

DLLS  840 

DLLS  850 

DLLS  860 

DLLS  920 

DLLS  930 

DLLS  97  0 

DLLS  980 

DLLS  990 

DLLS  1000 

DLLS  1010 

DLLS1020 

DLLS1030 

DLLS  1060 

DLLS  1070 

DLLS  1080 

DLL  SI  090 

DLL  SI  100 

DLLS1110 

DLLS  1120 

DLLS1130 

DLLS  1140 

DLLS  117  0 

DLLS  1180 

DLLS1190 

DLLS  12  00 

DLLS  12  10 

DLLS1240 

DLLS1270 

DLLS1280 

DLLS  1290 

DLLS  1320 

DLLS  13  30 

DLLS  1370 

DLLS  138  0 

DLLS  139  0 

DLLS1400 

DLLS  1410 

DLLS 1420 

DLLS  1450 

DLLS  3.460 

DLLS  147  0 

DLLS  1480 

DLLS1490 

DLLS1500 

DLLS  15  1.0 

DLLS  1520 

is 


16 


17 

18 
19 


2G 


21 


22 

2  3 


24 


25 
26 


27 
28 


I  I  =  I  +  T  0 

h  =  H  +  A(  I  )*A(  I  I  ) 
H=BETA*H 

00  16  1  =  1  ST, I  END 

II  =  1  +  10 

A( II  )=A(  II )-A(  I )*H 

II=IST+ID 

H=AUX< J )-A(  II  )*A(  I  I  ) 

AIIX(  J  )=H 

IF(H-PIV)  18,  18,  17 

PIV=H 

KPIV=J 

CONTINUE 

HO  21  J=.K,LM,"M 

M=O.DO 

1  E'ND=J+M-K 
1 1  =  i  S  T 

00  2  0  I=J,IEND 
H=H+A(  I  I  )*B(  I ) 

1  I  =  I  1  +  1 
H=BFTA*H 
11  =  1ST 

no  21  I=J,  I  END 

R( I )=B( I )-A( II )*H 

11=11+1 

IF^  =  0 

I=N 

LN=L*N 

PIV  =  1.00/AUX(  2*i\i) 

00  2  2  K=N,LN,N 
X(K )  =  PIV*R(  I  ) 

1  =  I+M 

IF(N-1)26,26,23  a 
JST=(N-1 )*M+N 

00  2  5  J  =  2,N 

JST=JST-M-1 

K=K!  +  M+1-J 

PIV=1.00/AIJX(  K  ) 

KST=K-M 

ID=IPIV(KST)-KST 

IST=2-J 

00  2  5  K  =  1,L 

H=B(KST) 

IST=IST+N 

IEN0=IST+J-2 

I I=JST 

00  24  1  =  1  ST,  I  END 

1  1  =  1 I+M 

H  =  H-A(  I  I  )*X(  I ) 

I = I ST-1 

I  I  =  I  +  10 

X(  I  )=X(  II  ) 

X{  I  I )=PIV*H 

KST=KST+M 

IST=N+1 

IEMD=0 

00  29  J=1,L 

1  END  =  I  END  +  M 
H=0.00 

IF( M-N)29,29, 2  7 
DO  28  1  =  1  ST,  I  END 
H=H+B( I )*B( I ) 


DLLS1530 
DLLS1540 


DLLS1550 
DLLS  156  0_ 
DLLS 1570 

DLLS1580 


DLLS1610 
DLLS1620 


DLLS1630 
DLLS1640 
DLLS1650 
DLLS1660 
DLLS1670 
0LLS1700 


200- 


DLLS1710 

DLLS1720 

0LLS1730 

DLLS1740 

DLLS  1750 

DL  L  S  1 7_6  0 

DLLS  1770 

DLLS1780 

DLLS  1790 

DLLS1800 

DLLS  18 10 

DLLS1860 

DLLS  18  70 

0LLS1880 

DLLS  1890 

DLLS1900 

DLLS  1910 

DLL  SI  92  0_ 

DLLS  1930 

DLLS194Q 

"DLLS  1950 

DLLS1960 

DLLS  1970 

D_LLS1980_ 

DLLS  1990 

DLLS2000 

DLLS 20 10 

0LLS2020 

DLLS2030 

DLLS2040 

DLLS2050 

0LLS2060 

DLLS207O 

DLLS2080 

0LLS2  09  0 

0LLS2100 

DLLS 211 6 

0LLS2120 

OLLS2130 

DLLS2140 

DLLS2180 

DLLS2190 

OLLS2200 

DLLS2220 
DLLS2230 

DLL S2 2 40 
OLLS2250 


X 


IST=IST+M 

29 

AUX( J)=H 

RETURN' 

32 

IER=K-1 

RETURN 

END 

DLLS2260 
DLLS22  70 
DLLS2280 
DLLS2390 

DLLS2400 
DLLS2410 
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